Echange_contact_PolyMAC_CDO.cpp

/****************************************************************************
* Copyright (c) 2024, 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 <Echange_contact_PolyMAC_CDO.h>
#include <Champ_front_calc.h>
#include <Probleme_base.h>
#include <Champ_Uniforme.h>
#include <Schema_Euler_Implicite.h>
#include <Milieu_base.h>
#include <Domaine_PolyMAC_CDO.h>
#include <Equation_base.h>
#include <Champ_Elem_PolyMAC_CDO.h>
#include <Op_Diff_PolyMAC_CDO_Elem.h>
 
Implemente_instanciable(Echange_contact_PolyMAC_CDO,"Paroi_Echange_contact_PolyMAC_CDO",Echange_externe_impose);
 
Sortie& Echange_contact_PolyMAC_CDO::printOn(Sortie& s ) const
{
  return s << que_suis_je() << "\n";
}
 
Entree& Echange_contact_PolyMAC_CDO::readOn(Entree& s )
{
  Nom nom_bord;
  Motcle nom_champ;
  s >> nom_autre_pb_ >> nom_bord >> nom_champ >> h_paroi;
  T_autre_pb_.typer("Champ_front_calc");
  Champ_front_calc& ch=ref_cast(Champ_front_calc, T_autre_pb());
  ch.creer(nom_autre_pb_, nom_bord, nom_champ);
  le_champ_front.typer("Ch_front_var_instationnaire_dep");
  T_ext().fixer_nb_comp(1);
  return s ;
}
 
void Echange_contact_PolyMAC_CDO::completer()
{
  Echange_externe_impose::completer();
  T_autre_pb().associer_fr_dis_base(T_ext().frontiere_dis());
  T_autre_pb().completer();
}
 
int Echange_contact_PolyMAC_CDO::initialiser(double temps)
{
  if (!Echange_externe_impose::initialiser(temps))
    return 0;
 
  Champ_front_calc& ch=ref_cast(Champ_front_calc, T_autre_pb());
  const Equation_base&    o_eqn  = ch.equation();
  const Front_VF& fvf = ref_cast(Front_VF, frontiere_dis()), o_fvf = ref_cast(Front_VF, ch.front_dis());
  const Domaine_PolyMAC_CDO& o_domaine = ref_cast(Domaine_PolyMAC_CDO, ch.domaine_dis());
  const IntTab& o_f_e = o_domaine.face_voisins(), &o_e_f = o_domaine.elem_faces();
  int i, j, k, l, e, f, N = ch.nb_comp(), o_n_f = o_e_f.dimension(1);
  //Nom nom_racc1=frontiere_dis().frontiere().le_nom();
 
  //arrays related to the Le Potier / Mahamane stabilization (see Op_Diff_PolyMAC_CDO_Elem)
  stab_ = ref_cast(Op_Diff_PolyMAC_CDO_Elem, o_eqn.operateur(0).l_op_base()).stab_;
 
  h_imp_.typer("Champ_front_fonc");
  h_imp_->fixer_nb_comp(N);
  h_imp_->valeurs().resize(0, N);
  fvf.frontiere().creer_tableau_faces(h_imp_->valeurs());
 
  ch.initialiser(temps,domaine_Cl_dis().equation().inconnue());
 
  monolithic = sub_type(Schema_Euler_Implicite, o_eqn.schema_temps()) ?
               ref_cast(Schema_Euler_Implicite, o_eqn.schema_temps()).resolution_monolithique(o_eqn.domaine_application()) : 0;
  if (!monolithic) return 1; //rest not needed
  o_domaine.init_m2(), o_domaine.init_virt_ef_map();
 
  /* src(i) = (proc, j) : source de l'item i de mdv_elem_faces */
  IntTab src(0, 2);
  MD_Vector_tools::creer_tableau_distribue(o_domaine.mdv_elems_faces, src);
  for (i = 0; i < src.dimension_tot(0); i++) src(i, 0) = Process::me(), src(i, 1) = i;
  src.echange_espace_virtuel();
 
  /* o_proc, o_item -> processor/item of the element, then of the other faces used in W2(f,.) for each face of the boundary */
  DoubleTrav o_proc(0, o_n_f), o_item(0, o_n_f), proc(0, o_n_f), l_item(0, o_n_f);
  o_domaine.creer_tableau_faces(o_proc), o_domaine.creer_tableau_faces(o_item);
  fvf.frontiere().creer_tableau_faces(proc), fvf.frontiere().creer_tableau_faces(l_item);
  int c_max = 1; //maximum number of items/coefficients to exchange per face
  for (i = 0; i < o_fvf.nb_faces(); i++)
    {
      f = o_fvf.num_face(i), e = o_f_e(f, 0);
      for (j = 0; j < o_n_f && o_e_f(e, j) != f; ) j++; // f = o_e_f(e, j)
      o_proc(f, 0) = src(e, 0), o_item(f, 0) = src(e, 1); //element
      for (k = o_domaine.w2i(o_domaine.m2d(e) + j) + 1, l = 1; k < o_domaine.w2i(o_domaine.m2d(e) + j + 1); k++, l++) //W2 hors diag
        {
          int idx = o_domaine.nb_elem_tot() + o_e_f(e, o_domaine.w2j(k));
          o_proc(f, l) = src(idx, 0), o_item(f, l) = src(idx, 1);
        }
      for (c_max = std::max(c_max, l); l < o_n_f; l++) o_proc(f, l) = o_item(f, l) = -1; //fill remaining entries with -1
    }
  c_max = Process::mp_max(c_max);
 
  //projection onto the local boundary
  if (o_fvf.frontiere().que_suis_je() == "Raccord_distant_homogene")
    o_fvf.frontiere().trace_face_distant(o_proc, proc), o_fvf.frontiere().trace_face_distant(o_item, l_item);
  else o_fvf.frontiere().trace_face_local(o_proc, proc), o_fvf.frontiere().trace_face_local(o_item, l_item);
 
  //filling
  item.resize(fvf.nb_faces(), c_max), item = -1;
  for (i = 0; i < fvf.nb_faces(); i++)
    for (j = 0; j < c_max && l_item(i, j) >= 0; j++)
      if (proc(i, j) == Process::me()) item(i, j) = (int)std::lrint(l_item(i, j));                     //local item (real)
      else
        {
          if (o_domaine.virt_ef_map.count({{ (int) proc(i, j), (int) l_item(i, j) }}))   //local item (virtual)
          {
            item(i, j) = o_domaine.virt_ef_map.at({{ (int) proc(i, j),  (int) l_item(i, j) }});
          }
          else Process::exit(Nom("Echange_contact_PolyMAC_CDO: missing item opposite face ") + Nom(fvf.num_face(i)) + " on " + fvf.le_nom() + " ! Have you used Decouper_multi?");
        }
 
  //coeff: 1 coefficient more than item -> the one for the face itself (placed at the beginning)
  coeff.resize(0, 1 + item.dimension(1), N), delta_int.resize(0, N, 2),
               fvf.frontiere().creer_tableau_faces(coeff);
 
  if (stab_)
    {
      delta.resize(0, item.dimension(1), N);
      fvf.frontiere().creer_tableau_faces(delta_int), fvf.frontiere().creer_tableau_faces(delta);
    }
  coeffs_a_jour_ = 0, delta_a_jour_ = (stab_ ? 1 : 0);
  return 1;
}
 
void Echange_contact_PolyMAC_CDO::update_coeffs()
{
  if (coeffs_a_jour_) return;
  //objects on the other side: equation, domain, unknown (prefix o_), boundary (for trace_face_distant)
  Champ_front_calc& ch=ref_cast(Champ_front_calc, T_autre_pb());
  const Domaine_PolyMAC_CDO& o_domaine = ref_cast(Domaine_PolyMAC_CDO, ch.domaine_dis());
  const Equation_base&    o_eqn  = ch.equation();
  const Op_Diff_PolyMAC_CDO_Elem& o_op_diff = ref_cast(Op_Diff_PolyMAC_CDO_Elem, o_eqn.operateur(0).l_op_base());
  const Champ_Elem_PolyMAC_CDO& o_inc  = ref_cast(Champ_Elem_PolyMAC_CDO, ch.inconnue());
  const Front_VF&         o_fvf  = ref_cast(Front_VF, ch.front_dis());
  const IntTab& e_f = o_domaine.elem_faces();
  const DoubleVect& fs = o_domaine.face_surfaces(), &ve = o_domaine.volumes();
 
  //"nu_faces" arrays used by the diffusion operators on each side
  o_op_diff.update_nu(), o_op_diff.update_delta_int();
  int ne_tot = o_domaine.nb_elem_tot(), N = ch.nb_comp();
 
  //arrays at faces on the other side, to be transmitted via o_fr.trace_face_{distant,local}: filled only for faces of o_fr
  o_domaine.init_m2();
  DoubleTrav o_Text, o_Himp, o_coeff, o_delta_int, nu_ef(e_f.dimension(1), N);
  if (monolithic) o_coeff.resize(o_domaine.nb_faces(), 1 + item.dimension(1), N), o_delta_int.resize(o_domaine.nb_faces(), N, 2);
  else o_Text.resize(o_domaine.nb_faces(), N), o_Himp.resize(o_domaine.nb_faces(), N);
  for (int i = 0; i < o_fvf.nb_faces(); i++)
    {
      int f = o_fvf.num_face(i), e = o_domaine.face_voisins(f, 0), fb, j, k, n, i_f = 0, idx;
      for (j = 0; j < e_f.dimension(1) && (fb = e_f(e, j)) >= 0; j++)
        if (fb == f) i_f = j; //index of face f within element e
      o_op_diff.remplir_nu_ef(e, nu_ef);
 
      /* build the outgoing heat flux, putting the T_f part in H_imp and the rest in T_ext */
      for (j = o_domaine.w2i(o_domaine.m2d(e) + i_f), idx = 0; j < o_domaine.w2i(o_domaine.m2d(e) + i_f + 1); j++, idx++)
        for (fb = e_f(e, o_domaine.w2j(j)), n = 0; n < N; n++)
          {
            double fac = fs(fb) / ve(e) * o_domaine.w2c(j) * nu_ef(o_domaine.w2j(j), n);
            if (monolithic)
              o_coeff(f, idx ? idx + 1 : 0, n) += fs(f) * fac, o_coeff(f, 1, n) -= fs(f) * fac; //cote face, cote element
            else
              {
                if (f == fb) o_Himp(f, n) += fac;
                else o_Text(f, n) -= fac * o_inc.valeurs()(ne_tot + fb, n);
                o_Text(f, n) += fac * o_inc.valeurs()(e, n);
              }
          }
      for (n = 0; stab_ && monolithic && n < N; n++)
        for (k = 0; k < 2; k++) o_delta_int(f, n, k) = o_op_diff.delta_f_int(f, n, k);
      for (n = 0; !monolithic && n < N; n++)
        if (o_Himp(f, n) > 1e-10) o_Text(f, n) /= o_Himp(f, n); //convert to actual T_ext
      for (n = 0; h_paroi < 1e9 && n < N; n++) //account for the wall thermal resistance
        {
          double fac = h_paroi / (h_paroi + (monolithic ? o_coeff(f, 0, n) / fs(f) : o_Himp(f, n)));
          if (monolithic)
            {
              if (stab_) o_delta_int(f, n, 0) *= fac;
              for (j = 0; j < o_coeff.dimension(1); j++) o_coeff(f, j, n) *= fac;
            }
          else o_Himp(f, n) *= fac;
        }
    }
 
  //transmission: either via Raccord_distant_homogene, or simple copy
  if (o_fvf.frontiere().que_suis_je() == "Raccord_distant_homogene")
    {
      if (monolithic) o_fvf.frontiere().trace_face_distant(o_coeff, coeff);
      else o_fvf.frontiere().trace_face_distant(o_Text, T_ext().valeurs()), o_fvf.frontiere().trace_face_distant(o_Himp, h_imp_->valeurs());
      if (monolithic && stab_) o_fvf.frontiere().trace_face_distant(o_delta_int, delta_int);
    }
  else
    {
      if (monolithic) o_fvf.frontiere().trace_face_local(o_coeff, coeff);
      else o_fvf.frontiere().trace_face_local(o_Text, T_ext().valeurs()), o_fvf.frontiere().trace_face_local(o_Himp, h_imp_->valeurs());
      if (monolithic && stab_) o_fvf.frontiere().trace_face_local(o_delta_int, delta_int);
    }
 
  if (monolithic) coeff.echange_espace_virtuel();
  else T_ext().valeurs().echange_espace_virtuel(), h_imp_->valeurs().echange_espace_virtuel();
  if (monolithic && stab_) delta_int.echange_espace_virtuel();
  coeffs_a_jour_ = 1;
}
 
void Echange_contact_PolyMAC_CDO::update_delta() const
{
  if (!monolithic || delta_a_jour_) return; //already done
  const Champ_front_calc& ch=ref_cast(Champ_front_calc, T_autre_pb());
  const Front_VF& o_fvf = ref_cast(Front_VF, ch.front_dis());
  const Domaine_PolyMAC_CDO& o_domaine = ref_cast(Domaine_PolyMAC_CDO, ch.domaine_dis());
  const Equation_base&    o_eqn  = ch.equation();
  const Op_Diff_PolyMAC_CDO_Elem& o_op_diff = ref_cast(Op_Diff_PolyMAC_CDO_Elem, o_eqn.operateur(0).l_op_base());
  const IntTab& e_f = o_domaine.elem_faces();
 
  //filling
  o_op_diff.update_delta();
  int i, j, n, N = ch.nb_comp();
  DoubleTrav o_delta(o_domaine.nb_faces(), item.dimension(1), N);
  for (i = 0; i < o_fvf.nb_faces(); i++)
    {
      int f = o_fvf.num_face(i), fb, e = o_domaine.face_voisins(f, 0), i_f = 0, idx;
      for (j = 0; j < e_f.dimension(1) && (fb = e_f(e, j)) >= 0; j++)
        if (fb == f) i_f = j; //index of face f within element e
      for (n = 0; n < N; n++) o_delta(f, 0, n) = o_op_diff.delta_e(e, n);
      for (j = o_domaine.w2i(o_domaine.m2d(e) + i_f) + 1, idx = 1; j < o_domaine.w2i(o_domaine.m2d(e) + i_f + 1); j++, idx++) //skip the first coefficient (diagonal)
        for (n = 0; n < N; n++) o_delta(f, idx, n) = o_op_diff.delta_f(e_f(e, o_domaine.w2j(j)), n);
    }
 
  //transmission
  if (o_fvf.frontiere().que_suis_je() == "Raccord_distant_homogene") o_fvf.frontiere().trace_face_distant(o_delta, delta);
  else o_fvf.frontiere().trace_face_local(o_delta, delta);
  delta.echange_espace_virtuel();
  delta_a_jour_ = 1;
}
 
void Echange_contact_PolyMAC_CDO::mettre_a_jour(double temps)
{
  Echange_externe_impose::mettre_a_jour(temps);
  coeffs_a_jour_ = 0, delta_a_jour_ = (stab_ ? 0 : 1);
}