R134a media - differentiation problems #1575

modelica-trac-importer · 2017-01-15T01:01:25Z

Reported by johan.windahl on 26 Sep 2014 16:25 UTC
Hello,
I noted that the implementation of Modelica.Media.R134a.R134a_ph.setState_phX is not using the same structure as in most of the other Modelica.Media setState_XXX implementations where the state is calculated in one line.

Implementing the expression in one line has the advantage that a tool (Dymola) can differentiate even that the function itself do not specify an derivative annotation (which is not possible due to the record contains both reals and non-reals), see example below:

model TestDifferentiation
  // This example will not work with R134a due to the code in setState_phX is not written in one row
  // but it will work with Medium=Modelica.Media.Water.WaterIF97_ph 
   package Medium=Modelica.Media.R134a.R134a_ph;
   Medium.ThermodynamicState state;
   Modelica.SIunits.Density d;
   Modelica.SIunits.Mass m;
   Modelica.SIunits.Pressure p(stateSelect=StateSelect.always);
   Modelica.SIunits.SpecificInternalEnergy u;
   Modelica.SIunits.SpecificEnthalpy h(stateSelect=StateSelect.always);
   parameter Modelica.SIunits.Volume V=1;

equation 
   m=V*d;
   state=Medium.setState_phX(p,h,Medium.reference_X);
   u=Medium.specificInternalEnergy(state);
   d=Medium.density(state);
   der(m)=time;
   der(m*u)=time;
end TestDifferentiation;

In general I see a problem with that the Modelica.Media structure relies on several undocumented special constructs that are needed by the tool for an efficient implementation.

Migrated-From: https://trac.modelica.org/Modelica/ticket/1575

modelica-trac-importer · 2017-01-15T01:01:34Z

Comment by wischhusen on 3 Dec 2014 11:53 UTC
There is a general problem behind:

There is no finite derivative for Integer phase contained in the state record which might be essentially the problem from the numerical point of view. In the water medium model phase is set to be always 0 (not known) which helps to solve the problem, but does not deliver the result the user would expect.

modelica-trac-importer · 2017-01-15T01:01:44Z

Comment by hubertus on 7 Dec 2015 14:00 UTC
I don't think that this is the issue,afaik this can be fixed by using the noDerivative annotation on the Integer phase. Also: I don't see how the phase can be known in any arbitrary function call, so what is wrong about the handling in the water properties?

modelica-trac-importer · 2017-01-15T01:02:15Z

Modified by hubertus on 7 Dec 2015 14:03 UTC

modelica-trac-importer · 2017-01-15T01:02:35Z

Comment by wischhusen on 23 Dec 2015 09:55 UTC
It requires some more features to enable inverse use of setState functions.

beutlich · 2019-10-05T20:55:27Z

Though the prerequisites have changed with the Modelica specification v3.4 allowing mixed
Real and non-Real record derivatives (MCP-0028), this issue with TestDifferentiation failing to simulate still is apparent in current MSL (both maint/3.2.3 and master) and current Modelica tools.

@hubertus65 @thorade @wischhusen @jwindahlModelon Do you think you can fix it for MSL v4.0.0?

beutlich · 2019-10-06T07:57:22Z

Here's what I tried so far.

Changing Modelica.Media.R134a.R134a_ph.setState_phX to

diff --git a/Modelica/Media/R134a.mo b/Modelica/Media/R134a.mo
index e1f2434c6..41778bde9 100644
--- a/Modelica/Media/R134a.mo
+++ b/Modelica/Media/R134a.mo
@@ -289,15 +289,17 @@ package R134a "R134a: Medium model for R134a"
       Modelica.SIunits.SpecificEnthalpy hv=dewEnthalpy(sat) "Vapor enthalpy";
 
     algorithm
-      state.p := p;
-      state.h := h;
-
-      state.phase := if ((state.h < hl) or (state.h > hv) or (state.p >
-        R134aData.data.FPCRIT)) then 1 else 2;
-
-      state.d := density_ph(p, h);
-      state.T := temperature_ph(p, h);
-      annotation (Documentation(info="<html>
+      state := ThermodynamicState(
+          d=density_ph(
+            p,
+            h),
+          T=temperature_ph(
+            p,
+            h),
+          phase=if ((h < hl) or (h > hv) or (p > R134aData.data.FPCRIT)) then 1 else 2,
+          h=h,
+          p=p);
+      annotation (Inline=true, Documentation(info="<html>
 <p>This function should be used by default in order to calculate the thermodynamic state record used as input by many functions.</p>
 <p>
 Example:

does neither help to simulate the TestDifferentiation model in Dymola 2020 or SimulationX 4.0.

However, changing state.phase to always be zero leads to successful simulating model in SimulationX 4.0, whereas Dymola 2020 still gives the same error as before.

wischhusen · 2019-10-21T12:09:42Z

I tried to define a dummy derivative function getPhase_ph_der (always returning 0) of getPhase_ph as a test. Dymola always warns about incompatible types for the outputs (which is Integer in both cases). Changing the output type of phase and phase_der to Real is accepted. Would it not make sense to symbolically substract Integer inputs from the list of derivatives (when differantiating functions) in any case?

wischhusen · 2019-10-23T06:59:25Z

These are the functions:

function getPhase_ph "Number of phases by pressure and specific enthalpy"
  extends Modelica.Icons.Function;
  input AbsolutePressure p "Pressure";
  input SpecificEnthalpy h "Specific enthalpy";

  output Integer phase "Number of phases";

protected 
  SaturationProperties sat(psat=p, Tsat=0)
    "Saturation temperature and pressure";
  Modelica.SIunits.SpecificEnthalpy hl=bubbleEnthalpy(sat)
    "Liquid enthalpy";
  Modelica.SIunits.SpecificEnthalpy hv=dewEnthalpy(sat) "Vapor enthalpy";

algorithm 
  //phase :=0;
  phase := if ((h < hl) or (h > hv) or (p > R134aData.data.FPCRIT)) then 1 else 2;

  annotation (derivative=getPhase_ph_der, Documentation(info="<html>
This function computes the number of phases for R134a depending on the inputs for absolute pressure and specific enthalpy. It makes use of cubic spline functions for liquid and vapor specific enthalpy.
</html>"));
end getPhase_ph;

function getPhase_ph_der
  "Derivative of number of phases by pressure and specific enthalpy"
  extends Modelica.Icons.Function;
  input AbsolutePressure p "Pressure";
  input SpecificEnthalpy h "Specific enthalpy";
  input Real p_der "Pressure derivative";
  input Real h_der "Specific enthalpy derivative";

  output Integer phase_der "Derivative of number of phases";

algorithm 
  phase_der := 0;

end getPhase_ph_der;

The error message received from Dymola 2020 is:

Check of Modelica.Media.R134a.R134a_ph.getPhase_ph:

 function Modelica.Media.R134a.R134a_ph.getPhase_ph specified derivative Modelica.Media.R134a.R134a_ph.getPhase_ph_der, but the outputs did not match.

Check of Modelica.Media.R134a.R134a_ph.getPhase_ph successful.

WARNINGS have been issued.

HansOlsson · 2019-10-23T08:22:16Z

I tried to define a dummy derivative function getPhase_ph_der (always returning 0) of getPhase_ph as a test. Dymola always warns about incompatible types for the outputs (which is Integer in both cases). Changing the output type of phase and phase_der to Real is accepted. Would it not make sense to symbolically substract Integer inputs from the list of derivatives (when differantiating functions) in any case?

The underlying problem are expressions that are differentiated as part of index reduction that also have discontinuities (in the original expression). For these cases an integer-expression can only change discontinuously and since it is the phase it seems likely that it actually changes.

So, having getPhase_ph_der does not seem good.

However, it may be that the original problem can be avoided - but it requires a bit more investigation, and e.g. using noDerivative does make sense since phase-changes are connected to having phX at phase-boundaries - you don't change phase away from phase-boundaries, right?

beutlich · 2019-10-23T09:25:31Z

but it requires a bit more investigation

Who can do it? @HansOlsson @wischhusen @hubertus65 @thorade ?

wischhusen · 2019-10-23T10:03:42Z

@HansOlsson For calculating the derivatives of the states we do not need the derivative of "phase". I did not find a better way to introduce the dummy function to make it clear.

And yes, the number of phases only change on the phase boundaries.

Well, the solution could be to introduce function setState_phX_der? And use noDerivative on state.phase, right?

wischhusen · 2019-10-23T10:39:50Z

Here is what I tried to do:

function setState_phX_der
  "Set state for pressure and specific enthalpy (X not used since single substance)"
  input AbsolutePressure p "Pressure";
  input SpecificEnthalpy h "Specific enthalpy";
  input MassFraction X[:]=reference_X "Mass fractions";
  input FixedPhase phase=0
    "2 for two-phase, 1 for one-phase, 0 if not known";
  input Real p_der;
  input Real h_der;
  output ThermodynamicStateDerivatives state_der "Thermodynamic state record";
algorithm 
  state_der.p_der := p_der;
  state_der.h_der := h_der;

  state_der.d_der := rho_ph_der(p, h, derivsOf_ph(
          p,
          h,
          getPhase_ph(p, h)), p_der, h_der);
  state_der.T_der := T_ph_der(p, h, derivsOf_ph(
          p,
          h,
          getPhase_ph(p, h)), p_der, h_der);
end setState_phX_der;

record ThermodynamicStateDerivatives
  Real p_der;
  Real h_der;
  Real d_der;
  Real T_der;
  // phase and X is skipped here!
end ThermodynamicStateDerivatives;

redeclare function extends setState_phX
  "Set state for pressure and specific enthalpy (X not used since single substance)"

protected 
  SaturationProperties sat(psat=p, Tsat=0)
    "Saturation temperature and pressure";
  Modelica.SIunits.SpecificEnthalpy hl=bubbleEnthalpy(sat)
    "Liquid enthalpy";
  Modelica.SIunits.SpecificEnthalpy hv=dewEnthalpy(sat) "Vapor enthalpy";

algorithm 
  state.p := p;
  state.h := h;

  state.phase := if ((state.h < hl) or (state.h > hv) or (state.p >
    R134aData.data.FPCRIT)) then 1 else 2;

  state.d := density_ph(p, h);
  state.T := temperature_ph(p, h);
  annotation (derivative(noDerivative=phase, noDerivative=X) = setState_phX_der);
end setState_phX;

The tool throws a warning that the outputs of both functions do not match.

HansOlsson · 2019-10-23T13:22:13Z

First a general observation: writing smoothOrder in the setState-function should ideally suffice. There should not be a need for the user to write derivatives if they can be auto-generated. That could fairly easily be supported (well, at least in updated tools), but it uncovers a worse problem.

The worse problem is that:
The phase-handling in the Media-libraries doesn't really work in Modelica, since state.phase is sort of a non-discrete integer variables and those shouldn't exist in Modelica (since 1.3?)

However, in many cases this is avoided in one of the following ways:

The phase is just 0 (meaning don't know), and thus doesn't change.
There is no explicit state-variable (and thus no state.phase) in the model and instead the state is directly sent to a function.
E.g., I did barely see any references to state.phase in R134a.

wischhusen · 2019-10-23T14:28:20Z

I think there is no essential reason, why we see phase in the state record. Actually, it is just an informational output. In contrast to the other thermodynamic states it can not be used to calculated other state variables from it. Thus it would be sufficient to place and calculated it in the BaseProperties model. But this change would be non-backward-compatible.

HansOlsson · 2019-10-30T12:48:48Z

I think there is no essential reason, why we see phase in the state record. Actually, it is just an informational output. In contrast to the other thermodynamic states it can not be used to calculated other state variables from it. Thus it would be sufficient to place and calculated it in the BaseProperties model. But this change would be non-backward-compatible.

Even if non-backward-compatible it may be the best solution. The question is how large impact it would have, i.e. whether users use phase from this record.

BTW: I noticed that phase seem different for this media compared to how I thought 2-phase water operated, and I would say that it is number of phases - not the phase.

wischhusen · 2019-11-01T10:56:31Z

Good question ... . I think, most users would use the functions and not the record. Thus, we need to convert these calls.

From the developers point of view we could substitute state.phase by evaluating input states like p, h, d, T in order to find out wether one is in one-phase or two-phase. I think that the idea was to get rid of these evaluations requiring some calls of saturation property functions (e.g. h_liq and h_vap) in the underlying function calls again and again. Perhaps, we will sacrifize computation speed when removing phase as an input at all.

Probably it would be better to have an additional input variable piped through from the top-level functions or models, right? In such a case we could use the noDerivative attribute and tools could do the index reduction. Would that be worth a try?

For a quick change, we could use the workaround of having a single line in the setState functions.

beutlich · 2020-03-02T19:43:26Z

Removing milestone.

wischhusen · 2020-03-10T09:16:19Z

The following model can only be translated if I change R134a.R134a_ph.setState_phX in the following way:

redeclare function extends setState_phX
"Set state for pressure and specific enthalpy (X not used since single substance)"
algorithm
state := ThermodynamicState(phase=0, p=p, h=h, d=density_ph(p, h), T=temperature_ph(p, h));
end setState_phX;

model R134aTest
package SI=Modelica.Units.SI;
extends Modelica.Icons.Example;
replaceable package Medium = Modelica.Media.R134a.R134a_ph "Medium model";
SI.SpecificEnthalpy h = 420e3;
SI.Pressure p = 1e5;
Medium.ThermodynamicState state = Medium.setState_phX(p, h);
SI.SpecificEnthalpy h_dew = Medium.dewEnthalpy(Medium.SaturationProperties(Tsat=Medium.saturationTemperature(state.p), psat=Medium.pressure(state)));
end ;

The change still resolves for the previous differentiation problem stated above. Nevertheless, it is not optimal.

dietmarw · 2020-03-10T09:29:46Z

Just since I see this. Please do not use this construct:

package SI=Modelica.Units.SI;

The proper way is:
import Modelica.Units.SI;

modelica-trac-importer assigned hubertus65 Jan 15, 2017

modelica-trac-importer added bug Critical/severe issue L: Media Issue addresses Modelica.Media labels Jan 15, 2017

modelica-trac-importer assigned wischhusen and unassigned hubertus65 Jan 15, 2017

modelica-trac-importer added this to the MSL3.2.2 milestone Jan 15, 2017

modelica-trac-importer modified the milestones: MSL_next-MINOR-version, MSL3.2.2 Jan 15, 2017

beutlich modified the milestones: MSL3.2.3, MSL_next-MAJOR-version Jun 24, 2018

slamer59 mentioned this issue Jul 6, 2018

Failed to differentiate the equation ExternalMedia.Media.CoolPropMedium.setState_ph_ modelica-3rdparty/ExternalMedia#13

Closed

beutlich mentioned this issue Oct 22, 2019

Modelica.Media.R134a problem #3163

Closed

beutlich added the P: high High priority issue label Oct 23, 2019

MartinOtter self-assigned this Jan 9, 2020

MartinOtter removed their assignment Jan 9, 2020

beutlich removed this from the MSL4.0.0 milestone Mar 2, 2020

dietmarw closed this as completed Mar 10, 2020

dietmarw reopened this Mar 10, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

R134a media - differentiation problems #1575

R134a media - differentiation problems #1575

modelica-trac-importer commented Jan 15, 2017 •

edited by beutlich

Loading

modelica-trac-importer commented Jan 15, 2017

modelica-trac-importer commented Jan 15, 2017

modelica-trac-importer commented Jan 15, 2017

modelica-trac-importer commented Jan 15, 2017

beutlich commented Oct 5, 2019

beutlich commented Oct 6, 2019

wischhusen commented Oct 21, 2019

wischhusen commented Oct 23, 2019

HansOlsson commented Oct 23, 2019

beutlich commented Oct 23, 2019

wischhusen commented Oct 23, 2019

wischhusen commented Oct 23, 2019 •

edited by beutlich

Loading

HansOlsson commented Oct 23, 2019

wischhusen commented Oct 23, 2019 •

edited

Loading

HansOlsson commented Oct 30, 2019

wischhusen commented Nov 1, 2019

beutlich commented Mar 2, 2020

wischhusen commented Mar 10, 2020

dietmarw commented Mar 10, 2020

R134a media - differentiation problems #1575

R134a media - differentiation problems #1575

Comments

modelica-trac-importer commented Jan 15, 2017 • edited by beutlich Loading

modelica-trac-importer commented Jan 15, 2017

modelica-trac-importer commented Jan 15, 2017

modelica-trac-importer commented Jan 15, 2017

modelica-trac-importer commented Jan 15, 2017

beutlich commented Oct 5, 2019

beutlich commented Oct 6, 2019

wischhusen commented Oct 21, 2019

wischhusen commented Oct 23, 2019

HansOlsson commented Oct 23, 2019

beutlich commented Oct 23, 2019

wischhusen commented Oct 23, 2019

wischhusen commented Oct 23, 2019 • edited by beutlich Loading

HansOlsson commented Oct 23, 2019

wischhusen commented Oct 23, 2019 • edited Loading

HansOlsson commented Oct 30, 2019

wischhusen commented Nov 1, 2019

beutlich commented Mar 2, 2020

wischhusen commented Mar 10, 2020

dietmarw commented Mar 10, 2020

modelica-trac-importer commented Jan 15, 2017 •

edited by beutlich

Loading

wischhusen commented Oct 23, 2019 •

edited by beutlich

Loading

wischhusen commented Oct 23, 2019 •

edited

Loading