Cahn_Hilliard_Convection.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 <Cahn_Hilliard_Convection.h>
#include <Probleme_base.h>
#include <Schema_Cahn_Hilliard_Navier_Stokes.h>
 
Implemente_instanciable(Cahn_Hilliard_Convection,"Cahn_Hilliard_Convection",Cahn_Hilliard);
 
Sortie& Cahn_Hilliard_Convection::printOn(Sortie& os) const
{
  return Cahn_Hilliard::printOn(os);
}
 
Entree& Cahn_Hilliard_Convection::readOn(Entree& is)
{
  return Cahn_Hilliard::readOn(is);
}
 
void Cahn_Hilliard_Convection::set_param(Param& param) const
{
  Cahn_Hilliard::set_param(param);
  // Spécifiques à l'équation de Cahn-Hilliard advective :
  param.ajouter_non_std("convection",(this), Param::REQUIRED);
  param.ajouter_condition("is_read_convection","The convection operator must be read, select negligeable type if you want to neglect it.");
}
 
int Cahn_Hilliard_Convection::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  if (mot=="convection")
    {
      Cerr << "Reading and typing of the convection operator : " << finl;
      const Champ_base& vitesse_transportante = vitesse_pour_transport();
      associer_vitesse(vitesse_transportante);
      terme_convectif_.associer_vitesse(vitesse_transportante);
      is >> terme_convectif_;
      return 1;
    }
  else
    return Cahn_Hilliard::lire_motcle_non_standard(mot,is);
}
 
const Champ_base& Cahn_Hilliard_Convection::vitesse_pour_transport() const
{
  return probleme().get_champ("vitesse");
}
 
/*! @brief Renvoie le nombre d'operateurs de l'equation: Pour Cahn-Hilliard c'est 3.
 *
 * @return (int) le nombre d'operateur de l'equation
 */
int Cahn_Hilliard_Convection::nombre_d_operateurs() const
{
  return 3;
}
 
int Cahn_Hilliard_Convection::nombre_d_operateurs_tot() const
{
  return nombre_d_operateurs();
}
 
/*! @brief Renvoie le i-eme operateur de l'equation:
 *
 * @param (int i) l'index de l'operateur a renvoyer
 * @return (Operateur&) l'operateur indexe par i
 */
const Operateur& Cahn_Hilliard_Convection::operateur(int i) const
{
  switch(i)
    {
    case 0:
      return terme_mobilite_;
    case 1:
      return terme_kappa_;
    case 2:
      return terme_convectif_;
    default :
      Cerr << "Error for Cahn_Hilliard_Convection::operateur(int i)" << finl;
      Cerr << "Cahn_Hilliard has " << nombre_d_operateurs() <<" operators "<<finl;
      Cerr << "and you are trying to access the " << i <<" th one."<< finl;
      exit();
    }
  // Pour les compilos!!
  return terme_mobilite_;
}
 
/*! @brief Renvoie le i-eme operateur de l'equation:
 * @param (int i) l'index de l'operateur a renvoyer
 * @return (Operateur&) l'operateur indexe par i
 */
Operateur& Cahn_Hilliard_Convection::operateur(int i)
{
  switch(i)
    {
    case 0:
      return terme_mobilite_;
    case 1:
      return terme_kappa_;
    case 2:
      return terme_convectif_;
    default :
      Cerr << "Error for Cahn_Hilliard_Convection::operateur(int i)" << finl;
      Cerr << "Cahn_Hilliard has " << nombre_d_operateurs() <<" operators "<<finl;
      Cerr << "and you are trying to access the " << i <<" th one."<< finl;
      exit();
    }
  // Pour les compilos!!
  return terme_mobilite_;
}
 
/*! @brief Returns the time derivative of the unknown I of the equation: dI/dt = M-1*(sum(operators(I) + sources))
 *
 *  In : derivee contains I (but immediatly set to 0)
 *  Out: derivee contains dI/dt
 *
 */
DoubleTab& Cahn_Hilliard_Convection::derivee_en_temps_inco(DoubleTab& derivee)
{
  derivee = 0.;
 
  DoubleTrav secmem(derivee);
  // secmem = sum(operators) + sources + equation specific terms
 
  if (implicite_==0)
    {
      // On cherche à construire : dc/dt = M-1 * (D mutilde^n)
      // avec mutilde^n = beta*dW/dc^n - Kc^n
      // où D est la matrice de l'opérateur lié à mu (terme_mobilite_)
      // et K est la matirce de l'opérateur lié à kappa (terme_kappa_)
 
      // On calcule d'abord mutilde^n = beta*dW/dc^n - Kc^n (fonction compute_mutilde)
      compute_mutilde();
      // Ici normalement les valeurs de mutilde (qui est l'inconnue de l'opérateur 0) sont mises à jour
      operateur(0).calculer(secmem); // secmem = D mutilde^n
      operateur(2).ajouter(secmem); // secmem += advection term
      les_sources.ajouter(secmem); // secmem += sources
      solveur_masse->appliquer(secmem); // M-1 * secmem
      derivee = secmem;
      derivee.echange_espace_virtuel();
 
    }
  else if (implicite_>0)
    {
      Cerr << "[Cahn_Hilliard_Convection] Warning: Be careful with implicit part in derivative computation." << finl;
      Process::exit();
    }
  else
    {
      Cerr << "Error in Cahn_Hilliard_Convection::derivee_en_temps_inco" << finl;
      Cerr << "implicite_ = " << implicite_ << " has not been initialized!" << finl;
      Cerr << "May be " << que_suis_je() << "::completer() method doesn't call Cahn_Hilliard::completer()" << finl;
      Process::exit();
    }
 
  return derivee;
}
 
DoubleTab Cahn_Hilliard_Convection::fonction_residu(const DoubleTab& c)
{
  DoubleTab residu = Cahn_Hilliard::fonction_residu(c);
  //
  const int nb_elem = domaine_dis().nb_elem();
  const int nb_param = nb_parametres_d_ordre();
  const double dt=schema_temps().pas_de_temps();
  const double theta=ref_cast(Schema_Cahn_Hilliard_Navier_Stokes, schema_temps()).coeff_theta();
  //
  DoubleTrav advection(concentration_->valeurs());
  terme_convectif_.calculer(c, advection);
  solv_masse().appliquer(advection);
  for (int elem = 0; elem < nb_elem; elem++)
    for (int p = 0; p < nb_param; p++)
      residu(elem, p) -= dt*theta*advection(elem,p);
  //
  residu.echange_espace_virtuel();
  return residu;
}