HS_INPUTS in tabular interpolation

[Marcoliveti97]

Hi everybody!
Do you know why the enthalpy-entropy pair of input is not supported for the tabular backend?
It one of the most common pair for many problems, thanks to the total enthalpy conservation, but unfortunately it seems to be unavailable. Do you think it could be soon implemented?

Thank you for your work

[ibell]

Because H,S, are always the most complicated inputs to handle. Neither input is an independent variable of the EOS and as a consequence, you need to iterate. That is even true in the case of the tabular interpolation. Someone could add the iteration, but it won't be me. Pull requests welcome.

[Marcoliveti97]

Hi Ian and thanks for your reply.
I could help implementing the routine, but my lack of proficiency in C++ would make it extremely inefficient.
If you know somebody among the developer interested in the topic I could work with him in this implementation.
Let me know!

[ibell]

Maybe @msaitta-mpr ?

[msaitta-mpr]

Is your idea to make a new h/s table, or perform a 1D solver on the existing p/h table?

[Marcoliveti97]

I don't really know how the tables are originally created, but I'd say that if the table is generated through p and h, a double interpolation should be enough to get a value from a given h and s. Am I incorrect?

[msaitta-mpr]

It should be possible to simply iterate on s, but there will be a reduction in speed. Often the reason that people use the tables is for the increased property lookup speed. I wasn't sure what you had in mind.

I can take a look into implementing an h/s input pair routine that uses such a routine.

[Marcoliveti97]

That's true, but I don't expect it to be slower than calling the EOS with h and s as inputs.
This would however allow consistency in codes where look-up tables are used, avoiding to call EOS (or PropsSI) whenever hs are used as inputs.

[Marcoliveti97]

Hello!
Is there any update on the issue? Do you think is it something that will be added to the code?

Thank you :)

[msaitta-mpr]

I have been playing around with the new code. However, I am still working to implement a solver for the saturation curves. Currently, they are tabulated by P/T, so a solver needs to be implemented to determine which phase (or two-phase) the state point is in. There are similar methods in the HEOS module, but I will need some more time to move one of them over and make it compatible.

[Marcoliveti97]

Excellent! Thank you for your time and dedication and for considering this request. Let me know if I can be of help in any way or if there is any news!

Cheers

[plevisse-naarea]

Hey everyone,
Have there been updates on this ? I'm developing a cycle code that uses a lot of H-S routines with the tabular interpolation for speed. I've tried looking into coding the iteration on s, but I'm not very familiar enough with what's under the hood to be very productive at it...

[ibell]

CoolProp is an open-source library supported by the community. Someone needs to go ahead and dig in!

…

On Fri, Aug 9, 2024 at 2:19 PM plevisse-naarea ***@***.***> wrote: Hey everyone, Have there been updates on this ? I'm developing a cycle code that uses a lot of H-S routines with the tabular interpolation for speed. I've tried looking into coding the iteration on s, but I'm not very familiar enough with what's under the hood to be very productive at it... — Reply to this email directly, view it on GitHub <#2318 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAOGC2YISFL3FGWY4KLDAXDZQUB37AVCNFSM6AAAAABMI5RRP2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDENZYGQ4TOMRZHE> . You are receiving this because you commented.Message ID: ***@***.***>

[plevisse-naarea]

Hi,
Took me some time yesterday but it seems to work not too badly :
` Inside the update function defined in TabularBackends.cpp

case HmolarSmolar_INPUTS :{

         std::size_t iL = 0, iV = 0;
           // A bit trickier than the HEOS case. Iterate on P
          // in order to find s

          _hmolar  = val1;
          _smolar  = val2;
          class Residual : public FuncWrapper1D
          {
          public:
              CoolProp::TabularBackend& HEOS;
              CoolPropDbl hmolar, smolar;
              PureFluidSaturationTableData& pure_saturation = HEOS.dataset->pure_saturation;
              PhaseEnvelopeData& phase_envelope = HEOS.dataset->phase_envelope;
              SinglePhaseGriddedTableData& single_phase_logph = HEOS.dataset->single_phase_logph;
              SinglePhaseGriddedTableData& single_phase_logpT = HEOS.dataset->single_phase_logpT;
              
              Residual(CoolProp::TabularBackend& HEOS,CoolPropDbl hmolar_spec, CoolPropDbl smolar_spec)
              : HEOS(HEOS), hmolar(hmolar_spec), smolar(smolar_spec){};
              double call(double P) {
                                      CoolPropDbl  smolar_calc;
                                      HEOS._p = P;
                                      HEOS._hmolar = hmolar;    
                                      HEOS._smolar = smolar;
                                      
                                      HEOS.using_single_phase_table = true;  // Use the table !
                                      std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max(), iclosest = 0;
                                      CoolPropDbl hL = 0, hV = 0;
                                      SimpleState closest_state;
                                      bool is_two_phase = false;
                                      // If phase is imposed, use it, but only if it's single phase:
                                      //   - Imposed two phase still needs to determine saturation limits by calling pure_saturation.is_inside().
                                      //   - There's no speed increase to be gained by imposing two phase.
                                      if ((HEOS.imposed_phase_index == iphase_not_imposed) || (HEOS.imposed_phase_index == iphase_twophase)) {
                                          if (HEOS.is_mixture) {
                                              is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, HEOS._p, iHmolar, HEOS._hmolar, iclosest, closest_state);
                                          } else {
                                              is_two_phase = pure_saturation.is_inside(iP, HEOS._p, iHmolar, HEOS._hmolar, iL, iV, hL, hV);
                                          }
                                      }
                                      // Phase determined or imposed, now interpolate results
                                      if (is_two_phase) {
                                          HEOS.using_single_phase_table = false;
                                          HEOS._Q = (static_cast<double>(HEOS._hmolar) - hL) / (hV - hL);
                                          if (!is_in_closed_range(0.0, 1.0, static_cast<double>(HEOS._Q))) {
                                              throw ValueError(
                                              format("vapor quality is not in (0,1) for hmolar: %g p: %g, hL: %g hV: %g ", static_cast<double>(HEOS._hmolar), HEOS._p, hL, hV));
                                          } else {
                                              HEOS.cached_saturation_iL = iL;
                                              HEOS.cached_saturation_iV = iV;
                                              HEOS._phase = iphase_twophase;
                                          }
                                          if (HEOS.is_mixture) {
                                          HEOS.phase_envelope_sat(phase_envelope, iSmolar, iP, HEOS._p); }
                                          else {
                                          smolar_calc=pure_saturation.evaluate(iSmolar, HEOS._p, HEOS._Q, HEOS.cached_saturation_iL, HEOS.cached_saturation_iV);
                                          }
                                      } else {
                                          HEOS.selected_table = SELECTED_PH_TABLE;
                                          // Find and cache the indices i, j
                                          HEOS.find_native_nearest_good_indices(single_phase_logph, HEOS.dataset->coeffs_ph, HEOS._hmolar, HEOS._p, HEOS.cached_single_phase_i,
                                                                          HEOS.cached_single_phase_j);
                                          // Recalculate the phase
                                          HEOS.recalculate_singlephase_phase();
                                          smolar_calc=HEOS.evaluate_single_phase_phmolar(iSmolar, HEOS.cached_single_phase_i, HEOS.cached_single_phase_j);
                                      }
                  double r = smolar_calc - smolar;
                  return r;
              }
          }  resid(*this,_hmolar, _smolar);
          std::string errstr;
          // Find minimum pressure
          bool good_Pmin = false;
          double Pmin = this->AS->p_triple();
          double rmin;
          do {
              try {
                  rmin = resid.call(Pmin);
                  good_Pmin = true;
              } catch (...) {
                  Pmin = Pmin*1.5;
              }
              if (Pmin > this->AS->pmax()) {
                  throw ValueError("Cannot find good Pmin");
              }
          } while (!good_Pmin);

          // Find maximum temperature
          bool good_Pmax = false;
          double Pmax = this->AS->pmax() ;  // Just a little above, so if we use Tmax as input, it should still work
          double rmax;
          do {
              try {
                  rmax = resid.call(Pmax);
                  good_Pmax = true;
              } catch (...) {
                  Pmax = Pmax*0.99;
              }
              if (Pmax < Pmin) {
                  throw ValueError("Cannot find good Pmax");
              }
          } while (!good_Pmax);
          if (rmin * rmax > 0 && std::abs(rmax) < std::abs(rmin)) {
              throw CoolProp::ValueError(format("HS inputs correspond to temperature above maximum pressure of EOS [%g Pa]", this->AS->pmax()));
          }
          Brent(resid, Pmin, Pmax, DBL_EPSILON, 1e-10, 100);
          break;
      }`

Two open questions on this :

• I iterated on P because it felt more intuitive to use the PH tables directly since we have the H-S inputs --> The HEOS backend iterates on T... I don't know what the risks associated with this are.

• I'm no Git expert so I don't quite know what the procedures for checking the validity of the code (or the respect of the proper C++ code standards) are, help on this would be appreciated :)

  Cheers,

[msaitta-mpr]

Iteration on P seems like the natural thing to me - it is what I had started working on but never finished. Iteration on T is done because HEOS is given in terms of RT as opposed to PH.

Guidelines for contributing are provided at:
https://github.com/CoolProp/CoolProp/blob/master/.github/CONTRIBUTING.md#your-first-code-contribution
https://github.com/CoolProp/CoolProp/wiki
Hopefully those can guide you through the initial git actions.