modelica-3rdparty / ExternalMedia

The ExternalMedia library provides a framework for interfacing external codes computing fluid properties to Modelica.Media-compatible component models.
53 stars 36 forks source link

Cannot compute arbitrary partial derivatives with function `partialDeriv_state` #88

Open ferrucci-franco opened 1 year ago

ferrucci-franco commented 1 year ago

Hi there,

I'm trying to use the function partialDeriv_state to compute an arbitrary partial derivative. For that, I have created the following test model:

model TestPartialDeriv
  // Medium package inside model:
  package Medium = ExternalMedia.Examples.CO2CoolProp;
  // Parameters:
  parameter Medium.AbsolutePressure p_start = 20e5 "Presure";
  parameter Medium.Temperature T_start = Modelica.Units.Conversions.from_degC(0) "Initial temperature";
  parameter Medium.Temperature T_end = Modelica.Units.Conversions.from_degC(80) "Final temperature";
  // Computed parameters:
  parameter Medium.SpecificEnthalpy h_start = Medium.specificEnthalpy_pT(p=p_start, T=T_start, phase=1);
  parameter Medium.SpecificEnthalpy h_end = Medium.specificEnthalpy_pT(p=p_start, T=T_end, phase=1);
  // Variables:
  Medium.ThermodynamicState state "thermodynamic state";
  Medium.SpecificEnthalpy h "Medium specific enthalpy";
  Medium.Temperature T;
  Medium.Density d;
  Medium.AbsolutePressure p;
  //Real my_der;
  equation
    state = Medium.setState_ph(p=p_start, h=h);
    h = h_start + time*(h_end - h_start);
    T = Medium.temperature_ph(p=p,h=h); // same as 'state.T' or 'Medium.temperature(state)'
    d = Medium.density_ph(p=p,h=h); // same as 'd = state.d' or 'Medium.density(state)'
    p = state.p;
    //my_der = Medium.partialDeriv_state(of="d", wrt="T", cst="p", state=state);
    annotation(experiment(StartTime = 0, StopTime = 1, Tolerance = 0.0001, Interval = 0.0002));
end TestPartialDeriv;

In this code, if I uncomment the variable my_der as well as the line my_der = Medium.partialDeriv_state(of="d", wrt="T", cst="p", state=state);, OpenModelica creates the corresponing EXE file but the simulation stops with the following error:

terminate called after throwing an instance of 'CoolProp::CoolPropError<(CoolProp::CoolPropBaseError::ErrCode)4>'
  what():  Your input name [drhodT|p] is not valid in get_parameter_index (names are case sensitive)

I tried to look inside the C++ files TestPartialDeriv.c and TestPartialDeriv_functions.c to track the error. It looks like there is a problem when calling the function TwoPhaseMedium_partialDeriv_state_C_impl from the ExternalMediaDLL file (maybe the argument is wrong?).

Thanks in advance,

Franco

ferrucci-franco commented 1 year ago

Info: I think the error I showed before is thrown by line 650 of coolpropsolver.cpp (function CoolPropSolver::makeDerivString):

long iOutput = CoolProp::get_parameter_index(derivTerm.c_str());

Going into CoolProp code, the exception shown in Modelica comes from CoolProp function get_parameter_index() (see here). The exception is thrown since function is_valid_parameter() (see is_valid_parameter) returns 0. Apparently, is_valid_parameter() cannot handle an input like dpdT|h.

In any case, in order to compute the partial derivative, coolpropsolver.cpp then uses the function double CoolProp::AbstractState::keyed_output(...) (see line 624 of coolpropsolver.cpp), which, I think, cannot handle partial derivatives parameters as inputs (see function here).

Would it be possible for coolpropsolver.cpp to implement TwoPhaseMedium_partialDeriv_state_C_impl() in the following way?:

  1. Inside CoolPropSolver::partialDeriv_state(), call Coolprop function CoolProp::get_parameter_index() for parameters of, wrt and cst to get the corresponding valid parameters (listed in CoolProp::parameters, such as iT, iP, iHmass, etc.).
  2. Call CoopProp function CoolProp::AbstractState::first_partial_deriv() (see here).

Thanks again,

Franco

casella commented 1 year ago

@ferrucci-franco I won't have time to look into this until September. If you want to experiment and find a working solution, please open a pull request. Thanks!