State functions do not seem to work properly

[PelDePinda]

Hi,

I am quite new with Modelica, but when I am trying to simulate the CoolProp models, Dymola gives an error about not being able to find a function.
This is the error message I get when trying the run the standard example ExternalMedia.Test.CoolProp.Pentane.Pentane_hs, delivered with the ExternalMedia package.

Error: The following error was detected at time: 0
ExternalMedia error: Inputs in Brent [143.470000,656.500000] do not bracket the root.  Function values are [6936.718242,139426.497618]
The stack of functions is:
setState_hs_Unique25
setState_hs_Unique25(h, 1500.0, 0)
Error: Integrator failed to start model.

The other example ExternalMedia.Test.CoolProp.Pentane.Pentane_hs_state does work however.

I get a similar error when making use of media described in Coolprop in another model (heat exchanger), such as ExternalMedia.Examples.CO2CoolProp or any other media described for that matter.

Failed to differentiate the equation
evaporatorGas.organicFlow.fluidState[2] = ExternalMedia.Media.CoolPropMedium.setState_phX_Unique35 (sinkPressure.p0, evaporatorGas.organicFlow.h[2], {1.0}, 0);

in order to reduce the DAE index.

Cannot find differentiation function:
ExternalMedia.Media.CoolPropMedium.setState_phX_Unique35(sinkPressure.p0, evaporatorGas.organicFlow.h[2], {1.0}, 0)
with respect to time

Dymola is currently unable handle index reduction involving non-real equations.

Failed to reduce the DAE index.

I suppose that Dymola is looking for this function that is nowhere defined (with the addition Unique35 or Unique25 that is).

[PStrm]

Hi,
I'm having a similar problem. Basically, calling a time derivative of a variable that is connected to setState_xx seem to be unsolvable cause of the DAE index reduction (in my example, the variable in question is "DER") . I've read that you should be able to provide the derivative by including the "derivative" annotation to the problematic function, but I haven't managed to accomplish that.

here is the model, it uses R290 from external media

model Compressor_test_small
 replaceable package Medium_ref = R290;

  //Medium_ref.ThermodynamicState State_in;
  Modelica.Units.SI.SpecificEnthalpy h_inlet;
  Modelica.Units.SI.SpecificEnthalpy h_ideal;
  Modelica.Units.SI.SpecificEnthalpy h_real;
  Modelica.Units.SI.AbsolutePressure p1;
  Modelica.Units.SI.AbsolutePressure p2;
  Real DER;
  Modelica.Units.SI.SpecificEntropy S;
equation 

  DER=der(h_ideal);

 // h_inlet = 400000;
  h_inlet = 400000+100*time;

 // p1 = 500000+100*time;
  p1 = 500000;

  p2 = 1500000;
  h_real = h_ideal * 1.2;

  //h_ideal = 600000;
  S = Medium_ref.specificEntropy_ph(p1,h_inlet);
  h_ideal = Medium_ref.specificEnthalpy_ps(p2,S);
 // h_ideal=Medium_ref.isentropicEnthalpy(p2,State_in);
 //State_in = Medium_ref.setState_ph(p1, h_inlet);

end Compressor_test_small;

since the function "specificEnthalpy" uses "setState_ps" I get an error

Failed to differentiate the equation
h_ideal = TESTY.R290.specificEnthalpy_ps_state_Unique7(
1500000.0, 
S, 
TESTY.R290.setState_ps_Unique9(1500000.0, S, 0));

in order to reduce the DAE index.

Cannot find differentiation function:
TESTY.R290.setState_ps_Unique9(1500000.0, S, 0)
with respect to time

Failed to reduce the DAE index.

defining the input variables h_inlet and p1 as constants leads to no problems, but is unlikely to happen in my simulations.

Thanks in advance for any answer
Pavel

[thorade]

Dymola and probably other tools can not calculate the derivatives of functions that use records as input or output.
So this will fail to be differentiated: density(setState_pT(101325, 300))
but this will work to get the time derivative: density_pT(101325, 300)

[PStrm]

This approach is what I am using - the problem is (for me) in the internals of the function - in my example, I do not use the setState_ps function - but the function "specificEnthalpy_ps" (and "isentropicEnthalpy") does use it:

h := specificEnthalpy_ps_state(p=p, s=s, state=setState_ps(p=p, s=s, phase=phase));

It's not like I am calling the derivative of the record - I am calling the derivative of just the enthalpy function basically ("Medium_ref.specificEnthalpy_ps") - but the solver itself calls the derivative of the state record

I have partially solved this issue with using a numerical kind of derivative for the enthalpy

when sample(2,1) then x = h_in; t = time; approx_der = (x1-x2) / (t-pre(t)); x1=x; x2=pre(x); Flag=1; end when;

but this creates issues with model stability and simulation time due to the sampling.

Thanks in advance for any answer,
Pavel

[thorade]

Quick link to the code in question for specificEnthalpy_ps function:

ExternalMedia/Modelica/ExternalMedia/Media/BaseClasses/ExternalTwoPhaseMedium.mo

Lines 771 to 784 in 8f5df11

	redeclare replaceable function specificEnthalpy_ps
	"Return specific enthalpy from p and s"
	extends Modelica.Icons.Function;
	input AbsolutePressure p "Pressure";
	input SpecificEntropy s "Specific entropy";
	input FixedPhase phase = 0
	"2 for two-phase, 1 for one-phase, 0 if not known";
	output SpecificEnthalpy h "specific enthalpy";
	algorithm
	h := specificEnthalpy_ps_state(p=p, s=s, state=setState_ps(p=p, s=s, phase=phase));
	annotation (
	Inline = true,
	inverse(s=specificEntropy_ph(p=p, h=h, phase=phase)));
	end specificEnthalpy_ps;

You want to achieve two things for good speed: derivatives and caching, and you have the annotations Inline or LateInline and derivative to do it. The caching is achieved by collecting all setState_ps function calls with identical input (you need inlining to move it to equation section).

The trick that is missing currently: specificEnthalpy_ps should not use Inline, it should use LateInline instead, and you need the derivative annotation. That should tell the tool to first find and use the derivative, and then afterwards do the inlining to get caching and good speed.

The above is at least true for Dymola and pure Modelica code, it might be different with OpenModelica and external code (maybe CoolPRop does some kind of caching??).
Also read this SO question&answer: https://stackoverflow.com/q/55019892/874701

[thorade]

There are variants to achieve the same, and looking at density_ph, it seem that ExternalMedia is using a different approach than what I descriebed:

ExternalMedia/Modelica/ExternalMedia/Media/BaseClasses/ExternalTwoPhaseMedium.mo

Lines 392 to 430 in 8f5df11

	redeclare replaceable function density_ph "Return density from p and h"
	extends Modelica.Icons.Function;
	input AbsolutePressure p "Pressure";
	input SpecificEnthalpy h "Specific enthalpy";
	input FixedPhase phase = 0
	"2 for two-phase, 1 for one-phase, 0 if not known";
	output Density d "Density";
	algorithm
	d := density_ph_state(p=p, h=h, state=setState_ph(p=p, h=h, phase=phase));
	annotation (Inline = true);
	end density_ph;
	function density_ph_state "returns density for given p and h"
	extends Modelica.Icons.Function;
	input AbsolutePressure p "Pressure";
	input SpecificEnthalpy h "Enthalpy";
	input ThermodynamicState state;
	output Density d "density";
	algorithm
	d := density(state);
	annotation (
	Inline=false,
	LateInline=true,
	derivative(noDerivative=state)=density_ph_der);
	end density_ph_state;
	replaceable function density_ph_der "Total derivative of density_ph"
	extends Modelica.Icons.Function;
	input AbsolutePressure p "Pressure";
	input SpecificEnthalpy h "Specific enthalpy";
	input ThermodynamicState state;
	input Real p_der "time derivative of pressure";
	input Real h_der "time derivative of specific enthalpy";
	output Real d_der "time derivative of density";
	algorithm
	d_der := p_der*density_derp_h(state=state)
	+ h_der*density_derh_p(state=state);
	annotation (Inline=true);
	end density_ph_der;

It is quite some effort to test whether caching works as expected, probably need to generate flat Modelica output, and also do some benchmarking. After that, the same approach can be used on all functions.

[PStrm]

I created a simple model using a derivative of the result of the density_ph and it seems that it works fine.

So basically, I create a function "specificEnthalpy_ps_der" calculating the total derivative (using partialDeriv_state probably?) and add the derivative annotation to the "xxx_state" function(s).

Thanks for your answers truly, I'm a newbie in this part of Modelica.

Pavel