Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cannot compute arbitrary partial derivatives with function partialDeriv_state #88

Open
ferrucci-franco opened this issue Jul 19, 2023 · 2 comments

Comments

@ferrucci-franco
Copy link

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
Copy link
Author

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
Copy link
Collaborator

casella commented Jul 20, 2023

@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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants