DP_Impose.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 <Schema_Temps_base.h>
#include <Probleme_base.h>
#include <Equation_base.h>
#include <DP_Impose.h>
#include <Param.h>
 
/*! @brief Reads the specifications of an imposed Delta P from an input stream.
 *
 * @param is an input stream
 * @return the modified input stream
 * @throws not a field
 */
 
// XD type_perte_charge_deriv objet_lecture type_perte_charge_deriv INHERITS_BRACE not_set
 
// XD type_perte_charge_dp type_perte_charge_deriv dp NO_BRACE DP field should have 3 components defining dp, dDP/dQ, Q0
// XD attr dp_field field_base dp_field REQ the parameters of the previous formula (DP = dp + dDP/dQ * (Q - Q0)):
// XD_CONT uniform_field 3 dp dDP/dQ Q0 where Q0 is a mass flow rate (kg/s).
 
// XD type_perte_charge_dp_regul type_perte_charge_deriv dp_regul BRACE Keyword used to regulate the DP value in order
// XD_CONT to match a target flow rate. Syntax : dp_regul { DP0 d deb d eps e }
// XD attr DP0 chaine DP0 REQ Reference pressure drop value (possibly time‑dependent)
// XD attr deb chaine deb REQ target flow rate in kg/s
// XD attr alpha chaine alpha REQ alpha string
 
 
// XD DP_Impose source_base DP_Impose NO_BRACE Source term to impose a pressure difference according to the formula : DP
// XD_CONT = dp + dDP/dQ * (Q - Q0)
// XD attr aco chaine(into=["{"]) aco REQ Opening curly bracket.
// XD attr dp_type type_perte_charge_deriv dp_type REQ mass flow rate (kg/s).
// XD attr surface chaine(into=["surface"]) surface REQ not_set
// XD attr bloc_surface bloc_lecture bloc_surface REQ Three syntaxes are possible for the surface definition block: NL2
// XD_CONT For VDF and VEF: { X|Y|Z = location subzone_name } NL2 Only for VEF: { Surface surface_name }. NL2 For
// XD_CONT PolyMAC_CDO { Surface surface_name Orientation champ_uniforme }.
// XD attr acof chaine(into=["}"]) acof REQ Closing curly bracket.
Entree& DP_Impose::lire_donnees(Entree& is)
{
  Motcle acc_ouverte("{");
  Motcle motlu;
  is >> motlu;
  if (motlu != acc_ouverte)
    {
      Cerr << "We expected the keyword" << acc_ouverte << " instead of " << motlu << finl;
      Process::exit();
    }
  is >> motlu;
  regul_ = 0;
  if (motlu == "dp")
    is >> DP_;
  else if (motlu == "dp_regul")
    {
      regul_ = 1;
      Nom alpha_str, deb_str, dp_str;
      Param param("dp_regul");
      param.ajouter("DP0", &dp_str, Param::REQUIRED);
      param.ajouter("deb", &deb_str, Param::REQUIRED);
      param.ajouter("alpha", &alpha_str, Param::REQUIRED);
      param.lire_avec_accolades(is);
      f_DP_.setNbVar(1), deb_cible_.setNbVar(1), alpha_.setNbVar(1);
      f_DP_.setString(dp_str), deb_cible_.setString(deb_str), alpha_.setString(alpha_str);
      f_DP_.addVar("t"), deb_cible_.addVar("t"), alpha_.addVar("t");
      f_DP_.parseString(), deb_cible_.parseString(), alpha_.parseString();
    }
  else
    {
      Cerr << "We expected the keyword dp or dp_regul instead of " << motlu << finl;
      Process::exit();
    }
  return is;
}
 
void DP_Impose::mettre_a_jour(double temps)
{
  if (regul_) return;
  DP_->mettre_a_jour(temps);
}
 
void DP_Impose::update_dp_regul(const Equation_base& eqn, double deb, DoubleVect& bilan)
{
  if (!regul_) return;
  double t = eqn.probleme().schema_temps().temps_courant(), dt = eqn.probleme().schema_temps().pas_de_temps();
  f_DP_.setVar(0, t), deb_cible_.setVar(0, t), alpha_.setVar(0, t);
  double deb_cible = deb_cible_.eval();
  if (std::abs(deb_cible) > 1e-10)
    {
      const double eps = alpha_.eval(), error = (deb - deb_cible) / deb_cible;
      fac_regul_ -= dt * eps * error;
    }
 
  //for the monitoring file: only on the master, because Source_base::imprimer() sums over all processes
  if (!Process::me()) bilan(0) = f_DP_.eval() * fac_regul_, bilan(1) = deb, bilan(2) = deb_cible;
}