Navier_Stokes_phase_field.cpp

/****************************************************************************
* Copyright (c) 2021, 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 <Convection_Diffusion_Phase_field.h>
#include <Source_Con_Phase_field_base.h>
#include <Navier_Stokes_phase_field.h>
#include <Source_Con_Phase_field.h>
#include <Operateur_Diff_base.h>
#include <Milieu_Phase_field.h>
#include <Champ_Fonc_Tabule.h>
#include <Assembleur_base.h>
#include <Probleme_base.h>
#include <Discret_Thyd.h>
#include <Perf_counters.h>
#include <Constituant.h>
#include <Domaine_VF.h>
#include <TRUSTTrav.h>
#include <Domaine.h>
#include <EChaine.h>
#include <Param.h>
#include <Debog.h>
 
Implemente_instanciable(Navier_Stokes_phase_field,"Navier_Stokes_phase_field",Navier_Stokes_std);
// XD navier_stokes_phase_field navier_stokes_standard navier_stokes_phase_field INHERITS_BRACE Navier Stokes equation
// XD_CONT for the Phase Field problem.
 
Sortie& Navier_Stokes_phase_field::printOn(Sortie& is) const
{
  return Navier_Stokes_std::printOn(is);
}
 
Entree& Navier_Stokes_phase_field::readOn(Entree& is)
{
  Navier_Stokes_std::readOn(is);
 
  //Creation du terme source de gravite pour discretisation VDF
  if (milieu().a_gravite())
    {
      Source t;
      Source& so = les_sources.add(t);
      Cerr << "Creation of the buoyancy source term for the Navier_Stokes_phase_field equation" << finl;
 
      Nom disc = discretisation().que_suis_je();
      if (disc != "VDF")
        {
          Cerr << "Only VDF discretization is allowed." << finl;
          Process::exit();
        }
 
      Nom type_so = "Source_Gravite_PF_VDF";
      so.typer_direct(type_so);
      so->associer_eqn(*this);
    }
  return is;
}
 
const Champ_Don_base& Navier_Stokes_phase_field::diffusivite_pour_transport() const
{
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, milieu());
  if (mil.get_boussi() == 0)
    {
      if (mil.get_diff_boussi() == 0)
        return mil.mu();
      else if (mil.get_diff_boussi() == 1)
        return le_fluide->viscosite_dynamique();
      else
        {
          Cerr << "If the Boussinesq hypothesis is not done," << finl;
          Cerr << "the keyword viscosite_dynamique_constante must be used to indicate" << finl;
          Cerr << "if the dynamical viscosity is constant or c dependent." << finl;
          Cerr << "This selection must be done before reading the diffusion term" << finl;
          Cerr << "of the Navier_Stokes_phase_field equation." << finl;
          Process::exit();
        }
    }
  else if (mil.get_boussi() == 1)
    return le_fluide->viscosite_cinematique();
  else
    {
      Cerr << "The use of Boussinesq hypothesis must be indicated" << finl;
      Cerr << "by the keyword approximation_de_boussinesq specification." << finl;
      Cerr << "This selection must be done before reading the diffusion term" << finl;
      Cerr << "of the Navier_Stokes_phase_field equation." << finl;
      Process::exit();
    }
  //Pour compilation
  return Navier_Stokes_phase_field::diffusivite_pour_transport();
}
 
/*! @brief Dicretise l'equation.
 *
 */
void Navier_Stokes_phase_field::discretiser()
{
  Navier_Stokes_std::discretiser();
 
  // Selection de l'assembleur en pression specifique au Phase_Field (P_VDF de BM)
  Nom type = "";
  type = "Assembleur_P_VDF_Phase_Field";
  Cerr << "** Pressure assembling tool : " << type << " **" << finl;
  assembleur_pression_.typer(type);
  assembleur_pression_->associer_domaine_dis_base(domaine_dis());
  assembleur_pression_->set_resoudre_increment_pression(1);
}
 
double Navier_Stokes_phase_field::calculer_pas_de_temps() const
{
  // RLT: plutot que de multiplier dt(diff) par rho ici comme c'est fait dans TrioCFD 1.8.0 (cf. coeff ci-dessous)
  // (en considerant qu'il a ete evalue dans l'operateur de diffusion (Op_Diff_VDF_Face_base) avec la viscosite dynamique - ce qui n'est vrai que pour boussi_==0)
  // on associe le champ masse volumique a l'operateur de diffusion dans le cas boussi_==0 pour que le calcul de son dt soit correct directement
  double coeff = mp_min_vect(ref_cast(Milieu_Phase_field, milieu()).rho().valeurs());
 
  double dt = 0;
  double dt_op;
 
  Cerr << "" << finl;
  double tps_courant = schema_temps().temps_courant();
  double num_dt = schema_temps().nb_pas_dt();
  Cerr << "                           *--*" << finl;
  Cerr << "" << finl;
  Cerr << "Calculation at time : " << tps_courant << " s" << finl;
  Cerr << "Time step number : " << num_dt << finl;
  Cerr << "" << finl;
 
  dt = 1. / terme_convectif.calculer_pas_de_temps();
  if (le_schema_en_temps->limpr())
    {
      Cout << "NS_PF : pas de temps de convection : " << 1. / dt << finl;
    }
  dt_op = terme_diffusif.calculer_pas_de_temps() * coeff;
  if (le_schema_en_temps->limpr())
    {
      Cout << "NS_PF : pas de temps de diffusion : " << dt_op << finl;
    }
  dt = dt + 1. / dt_op;
  dt = std::min(1. / dt, le_schema_en_temps->pas_temps_max());
  if (le_schema_en_temps->limpr())
    {
      Cout << "NS_PF : pas de temps de calcul : " << dt << finl;
    }
 
  // Calcul du temps effectue en simulation, ainsi que de l'estimation du temps restant
  double t_total;
  double t_remaining;
  double tps_init = schema_temps().temps_init();
  double tps_max = schema_temps().temps_max();
 
  double t_passe=statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time);
  if (tps_courant!=tps_init)
    {
      t_total = t_passe * ((tps_max - tps_init) / (tps_courant - tps_init));
      t_remaining = t_total - t_passe;
      Cerr << "" << finl;
      Cerr << "**************** ELAPSED TIME FOR THE SIMULATION *******************" << finl;
      Cerr << "Elapsed time of calculation until current time : " << floor(t_passe) << " s" << finl;
      Cerr << "Estimated time needed for the remaining part of the calculation : " << floor(t_remaining) << " s" << finl;
      Cerr << "                            --" << finl;
      Cerr << "Estimated time for the complete calculation : " << floor(t_total) << " s" << finl;
      Cerr << "********************************************************************" << finl;
      Cerr << "" << finl;
    }
  return dt;
}
 
/*! @brief Interpole la masse volumique aux faces (pour l'assembleur pression en particulier)
 *
 */
void Navier_Stokes_phase_field::rho_aux_faces(const DoubleTab& tab_rho, Champ_Don_base& rho_face_P)
{
  const Domaine_VF& domaineVF = ref_cast(Domaine_VF, domaine_dis());
  const IntTab& face_voisins = domaineVF.face_voisins();
  const DoubleVect& volumes = domaineVF.volumes();
  int nbfaces = domaineVF.nb_faces();
  const int nbfaces_bord = domaineVF.premiere_face_int();
  int el0, el1;
  double vol0, vol1;
 
  // Calcul de rho_face_P aux faces pour l'assemblage avec l'assembleur pression PF (cf AssembleurPVDF de BM)
 
  for (int fac = 0; fac < nbfaces_bord; fac++)
    {
      el0 = face_voisins(fac, 0);
      if (el0 != -1)
        {
          rho_face_P.valeurs()(fac) = tab_rho(el0);
        }
      else
        {
          el1 = face_voisins(fac, 1);
          rho_face_P.valeurs()(fac) = tab_rho(el1);
        }
    }
 
  for (int fac = nbfaces_bord; fac < nbfaces; fac++)
    {
      el0 = face_voisins(fac, 0);
      el1 = face_voisins(fac, 1);
      vol0 = volumes(el0);
      vol1 = volumes(el1);
      rho_face_P.valeurs()(fac) = (vol0 * tab_rho(el0) + vol1 * tab_rho(el1)) / (vol0 + vol1);
    }
}
 
/*! @brief Dicretise l'equation.
 *
 */
int Navier_Stokes_phase_field::preparer_calcul()
{
  Cerr << "** Navier_Stokes_phase_field::preparer_calcul **" << finl;
  //   Cerr << " - Copie via le Front_Tracking et modifie pour le Phase_field -" << finl;
  Equation_base::preparer_calcul();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, milieu());
  if (mil.get_boussi() == 0)
    {
      OWN_PTR(Champ_Don_base) rho_face_P;
      rho_face_P.typer("Champ_Fonc_Face");
      rho_face_P->valeurs() = inconnue().valeurs();
      const DoubleTab& tab_rho = mil.rho().valeurs();
      rho_aux_faces(tab_rho, rho_face_P);
      assembleur_pression_->assembler_rho_variable(matrice_pression_, rho_face_P.valeur());
    }
  else
    assembleur_pression_->assembler(matrice_pression_);
 
  la_pression->changer_temps(schema_temps().temps_courant());
  la_pression->valeurs().echange_espace_virtuel();
 
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, gradP av", gradient_P->valeurs());
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, la_pression av", la_pression->valeurs());
  //Cerr << "Navier_Stokes_phase_field::preparer_calcul() avant grad" << finl;
  gradient.calculer(la_pression->valeurs(), gradient_P->valeurs());
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, gradP ap", gradient_P->valeurs());
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, la_pression ap", la_pression->valeurs());
 
  gradient_P->changer_temps(schema_temps().temps_courant());
  divergence_U->valeurs() = 0.;
  divergence_U->changer_temps(schema_temps().temps_courant());
  Cerr << "Projection of initial and boundaries conditions " << finl;
 
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, gradP av projeter", gradient_P->valeurs());
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, la_pression av projeter", la_pression->valeurs());
  if (projection_a_faire())
    projeter();
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, gradP ap projeter", gradient_P->valeurs());
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, la_pression ap projeter", la_pression->valeurs());
 
  Navier_Stokes_std::calculer_la_pression_en_pa();
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, vitesse", inconnue());
  Debog::verifier("Navier_Stokes_phase_field::preparer_calcul, pression", la_pression);
 
  if (le_traitement_particulier)
    le_traitement_particulier->preparer_calcul_particulier();
 
  return 1;
}
 
void Navier_Stokes_phase_field::mettre_a_jour(double temps)
{
  Navier_Stokes_std::mettre_a_jour(temps);
 
  Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, milieu());
 
  mil.update_rho_mu(temps);
 
  if (mil.get_boussi() == 0)
    {
      OWN_PTR(Champ_Don_base) rho_face_P;
      rho_face_P.typer("Champ_Fonc_Face");
      rho_face_P->valeurs() = inconnue().valeurs();
      const DoubleTab& tab_rho = mil.rho().valeurs();
      rho_aux_faces(tab_rho, rho_face_P);
      assembleur_pression_->assembler_rho_variable(matrice_pression_, rho_face_P.valeur());
    }
  else
    assembleur_pression_->assembler(matrice_pression_);
}
 
DoubleTab& Navier_Stokes_phase_field::derivee_en_temps_inco(DoubleTab& vpoint)
{
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, milieu());
  if (mil.get_boussi() == 0)
    {
 
      //   Cerr << "** Navier_Stokes_phase_field::derivee_en_temps_inco **" << finl;
      //   Cerr << " - Copie via Navier_Stokes_Front_Tracking et modifie pour le Phase_field - " << finl;
      DoubleTrav secmem(la_pression->valeurs());
      DoubleTrav gradP(la_vitesse->valeurs());
      double dt = le_schema_en_temps->pas_de_temps();
 
      //   const Probleme_base& pb=ref_cast(Probleme_base, probleme());
      const Domaine_VF& zvf = ref_cast(Domaine_VF, domaine_dis());
      const DoubleTab& tab_rho = mil.rho().valeurs();
      const IntTab& face_voisins = zvf.face_voisins();
      const int nbfaces_bord = zvf.premiere_face_int();
      const int nbfaces = zvf.nb_faces();
      //   const int nbelems = zvf.nb_elem_tot();
 
      // Calcul de rho sur les faces
      DoubleTab rho_face(nbfaces);
      int face, elem1, elem2, k;
      double rho_moyen;
      for (face = 0; face < nbfaces_bord; face++)
        {
          elem1 = face_voisins(face, 0);
          if (elem1 != -1)
            rho_face[face] = tab_rho[elem1];
          else
            {
              elem2 = face_voisins(face, 1);
              rho_face[face] = tab_rho[elem2];
            }
        }
      for (face = nbfaces_bord; face < nbfaces; face++)
        {
          // On generalise le calcule de la masse volumique
          elem1 = face_voisins(face, 0);
          elem2 = face_voisins(face, 1);
          rho_moyen = (zvf.volumes(elem1) * tab_rho[elem1] + zvf.volumes(elem2) * tab_rho[elem2]);
          rho_moyen /= zvf.volumes(elem1) + zvf.volumes(elem2);
          rho_face[face] = rho_moyen;
        }
 
      // 1/dt B Un :;
      divergence.calculer(la_vitesse->valeurs(), secmem);
      secmem /= dt;
 
      vpoint = 0.;
 
      // Ajout de Div(tau_d), terme diffusif
      DoubleTrav diff(vpoint);
      terme_diffusif.ajouter(la_vitesse->valeurs(), diff);
 
      if (mil.get_diff_boussi() == 1)
        {
          // on multiplie par rho :
          if (diff.line_size() == 1)
            {
              for (face = 0; face < nbfaces; face++)
                {
                  diff(face, 0) *= rho_face[face];
                  diff(face, 0) /= mil.rho0();
                }
            }
          else
            {
              double drho;
              for (face = 0; face < nbfaces; face++)
                {
                  drho = rho_face[face];
                  for (k = 0; k < dimension; k++)
                    {
                      diff(face, k) *= drho;
                      diff(face, k) /= mil.rho0();
                    }
                }
            }
          // ainsi on a : rho/rho0 div(tau_d), et quand on va diviser par rho ensuite, ne restera que div(tau_d)/rho0
        }
 
      vpoint += diff;
 
      // Modif pour convection :
      DoubleTrav conv(vpoint);
      conv = 0.;
      terme_convectif.ajouter(la_vitesse->valeurs(), conv);
 
      // on multiplie par rho :
      if (conv.line_size() == 1)
        {
          for (face = 0; face < nbfaces; face++)
            conv(face, 0) *= rho_face[face];
        }
      else
        {
          double drho;
          for (face = 0; face < nbfaces; face++)
            {
              drho = rho_face[face];
              for (k = 0; k < dimension; k++)
                conv(face, k) *= drho;
            }
        }
      vpoint += conv;
      // Ajout de rho u Grad(u), terme convectif
 
      // on divise vpoint par rho :
      if (vpoint.line_size() == 1)
        {
          for (face = 0; face < nbfaces; face++)
            vpoint(face, 0) /= rho_face[face];
        }
      else
        {
          double drho;
          for (face = 0; face < nbfaces; face++)
            {
              drho = rho_face[face];
              for (k = 0; k < dimension; k++)
                vpoint(face, k) /= drho;
            }
        }
 
      les_sources.ajouter(vpoint);
      vpoint.echange_espace_virtuel();
 
      solveur_masse->appliquer(vpoint);
      if (calculate_time_derivative())
        derivee_en_temps().valeurs() = vpoint;
      schema_temps().modifier_second_membre((*this), vpoint);
 
      // Appliquer le solveur masse => diviser par le volume
      // A la base TRUST gere des variables * volume. Donc ici on retombe sur les VRAIES variables.
 
      divergence.ajouter(vpoint, secmem);
      // On calcule Div(vpoint)/dt i-e Div(u_n+1)/dt et on l'ajoute a secmem qui comprend deja Div(u_n)/dt
 
      secmem *= -1;
      secmem.echange_espace_virtuel();
 
      // Resolution en pression :
      //Assembleur_P_VDF_FT& assembleur_pression=ref_cast(Assembleur_P_VDF_FT, assembleur_pression_.valeur());
      //assembleur_pression.modifier_matrice_pression(matrice_pression_);
 
      const Domaine_VF& domaine_VF = ref_cast(Domaine_VF, le_dom_dis.valeur());
      const DoubleVect& volumes = domaine_VF.volumes();
      int el0, el1;
      double vol0, vol1;
 
      // Calcul de rho_face_P aux faces pour l'assemblage avec l'assembleur pression PF (cf AssembleurPVDF de BM)
      //Champ_Fonc_Face rho_face_P;
      OWN_PTR(Champ_Don_base) rho_face_P;
      rho_face_P.typer("Champ_Fonc_Face");
      rho_face_P->valeurs() = inconnue().valeurs();
 
      for (int fac = 0; fac < nbfaces_bord; fac++)
        {
          el0 = face_voisins(fac, 0);
          if (el0 != -1)
            rho_face_P->valeurs()(fac) = tab_rho(el0);
          else
            {
              el1 = face_voisins(fac, 1);
              rho_face_P->valeurs()(fac) = tab_rho(el1);
            }
        }
 
      for (int fac = nbfaces_bord; fac < nbfaces; fac++)
        {
          el0 = face_voisins(fac, 0);
          el1 = face_voisins(fac, 1);
          vol0 = volumes(el0);
          vol1 = volumes(el1);
          rho_face_P->valeurs()(fac) = (vol0 * tab_rho(el0) + vol1 * tab_rho(el1)) / (vol0 + vol1);
        }
 
      assembleur_pression_->assembler_rho_variable(matrice_pression_, rho_face_P.valeur());
      assembleur_pression_->modifier_secmem(secmem);
 
      solveur_pression_->reinit();
      secmem.echange_espace_virtuel();
      solveur_pression_.resoudre_systeme(matrice_pression_.valeur(), secmem, la_pression->valeurs());
      assembleur_pression_->modifier_solution(la_pression->valeurs());
 
      la_pression->valeurs().echange_espace_virtuel();
 
      // M-1 Bt P
      gradient.calculer(la_pression->valeurs(), gradP);
      gradP.echange_espace_virtuel();
 
      solveur_masse->appliquer(gradP);
      // Idem qu'avant
 
      // Pression a diviser par rho :
      if (gradP.line_size() == 1)
        {
          for (face = 0; face < nbfaces; face++)
            gradP(face, 0) /= rho_face[face];
        }
      else
        {
          double drho;
          for (face = 0; face < nbfaces; face++)
            {
              drho = rho_face[face];
              for (k = 0; k < dimension; k++)
                gradP(face, k) /= drho;
            }
        }
 
      // Correction en pression
      vpoint -= gradP;
      vpoint.echange_espace_virtuel();
    }
 
  else
    Navier_Stokes_std::derivee_en_temps_inco(vpoint);
 
  return vpoint;
}