print and assert not executed at the right time in FMI-ME

Open casella opened this issue 1 year ago • 1 comments

Please consider the following test model

model Test "Model that checks the correct implementation of the 1st order derivative of the flow function"
  extends Modelica.Icons.Example;
  constant Real gain = 2 "Gain for computing the mass flow rate";
  parameter Real k = 0.35 "Flow coefficient";
  parameter Modelica.Units.SI.MassFlowRate m_flow_turbulent = 0.36 "Mass flow rate where transition to turbulent flow occurs";
  Modelica.Units.SI.MassFlowRate m_flow "Mass flow rate";
  Modelica.Units.SI.MassFlowRate m_flow_comp "Comparison value for m_flow";
  Modelica.Units.SI.PressureDifference dp "Pressure drop";
  Modelica.Units.SI.MassFlowRate err "Integration error";
initial equation
  m_flow = m_flow_comp;
equation
  dp = time*gain;
  m_flow = Buildings.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp(dp = dp, k = k, m_flow_turbulent = m_flow_turbulent);
  der(m_flow) = der(m_flow_comp);
  err = m_flow - m_flow_comp;
  //assert(abs(err) < 1E-3, "Error in implementation.");
  Modelica.Utilities.Streams.print("m_flow = "+String(m_flow, ".10f") + ", m_flow_comp = " + String(m_flow_comp, ".10f")+", time = " + String(time));
  annotation(
    experiment(StartTime = -2, StopTime = 2, Interval = 0.1, Tolerance = 1e-8),
    uses(Buildings(version = "11.0.1")));
end Test;

which uses the latest version on the master branch of Buildings from https://github.com/lbl-srg/modelica-buildings/. The model checks the correctness of the derivative of a function by comparing the function's output with the integral of its derivative, using a very tight tolerance 1e-8. The model is a slight variation of Buildings.Fluid.BaseClasses.FlowModels.Validation.BasicFlowFunction_dp_DerivativeCheck, without the assert.

Running this with OpenModelica I get the following output:

m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2 
The initialization finished successfully without homotopy method.
m_flow = -0.6062177826, m_flow_comp = -0.6062177009, time = -1.5 
m_flow = -0.4949747468, m_flow_comp = -0.4949745789, time = -1 
m_flow = -0.3499471450, m_flow_comp = -0.3499468676, time = -0.5 
m_flow = 0.0000000000, m_flow_comp = 0.0000003584, time = 0 
m_flow = 0.3499471450, m_flow_comp = 0.3499475075, time = 0.5 
m_flow = 0.4949747468, m_flow_comp = 0.4949751075, time = 1 
m_flow = 0.6062177826, m_flow_comp = 0.6062181302, time = 1.5 
m_flow = 0.7000000000, m_flow_comp = 0.7000003459, time = 2 
m_flow = 0.7000000000, m_flow_comp = 0.7000003459, time = 2

As expected, the result of the integration m_flow_comp is off by a small quantity, about 5e-7 relative error maximum.

If I export this model as an FMU with OpenModelica and run this Python script

from fmpy import *
fmu = './Test.fmu'
dump(fmu)
result = simulate_fmu(fmu,fmi_type='ModelExchange',output_interval = 0.5)

I get the following output:

LOG_STDOUT        | info    | m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2
LOG_STDOUT        | info    | m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2
LOG_STDOUT        | info    | m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2
LOG_STDOUT        | info    | m_flow = -0.5989604513, m_flow_comp = -0.6062176076, time = -1.5
LOG_STDOUT        | info    | m_flow = -0.6062177826, m_flow_comp = -0.6062176076, time = -1.5
LOG_STDOUT        | info    | m_flow = -0.4887926395, m_flow_comp = -0.4949744062, time = -1
LOG_STDOUT        | info    | m_flow = -0.4949747468, m_flow_comp = -0.4949744062, time = -1
LOG_STDOUT        | info    | m_flow = -0.3492669144, m_flow_comp = -0.3499465272, time = -0.5
LOG_STDOUT        | info    | m_flow = -0.3499471450, m_flow_comp = -0.3499465272, time = -0.5
LOG_STDOUT        | info    | m_flow = 0.0200833821, m_flow_comp = 0.0000006512, time = 0
LOG_STDOUT        | info    | m_flow = 0.0000000000, m_flow_comp = 0.0000006512, time = 0
LOG_STDOUT        | info    | m_flow = 0.3534517110, m_flow_comp = 0.3499478013, time = 0.5
LOG_STDOUT        | info    | m_flow = 0.3499471450, m_flow_comp = 0.3499478013, time = 0.5
LOG_STDOUT        | info    | m_flow = 0.4956841774, m_flow_comp = 0.4949754538, time = 1
LOG_STDOUT        | info    | m_flow = 0.4949747468, m_flow_comp = 0.4949754538, time = 1
LOG_STDOUT        | info    | m_flow = 0.6092723340, m_flow_comp = 0.6062183965, time = 1.5
LOG_STDOUT        | info    | m_flow = 0.6062177826, m_flow_comp = 0.6062183965, time = 1.5
LOG_STDOUT        | info    | m_flow = 0.7075988253, m_flow_comp = 0.7000005921, time = 2
LOG_STDOUT        | info    | m_flow = 0.7000000000, m_flow_comp = 0.7000005921, time = 2

This is a bit weird: the print statement is executed twice per reporting interval, and only the second instance shows the correct m_flow_comp off by a small quantity, while the first one has a much larger error, maybe it's the result at the previous variable step-size step or something.

What is more troublesome is that if uncomment the assert statement I get

LOG_STDOUT        | info    | m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2
LOG_STDOUT        | info    | m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2
LOG_STDOUT        | info    | m_flow = -0.7000000000, m_flow_comp = -0.7000000000, time = -2
LOG_STDOUT        | info    | m_flow = -0.5989604513, m_flow_comp = -0.6062176076, time = -1.5
LOG_ASSERT        | error   | [/home/casella/Modelicon/Test.mo:17:3-17:54:writable]
|                 | |       | The following assertion has been violated at time -1.500000
|                 | |       | ((abs(err) < 0.001)) --> "Error in implementation."

which means that the wrong value is actually used to test the assert, which obviously fails. This causes the failure of 13 test cases from the Buildings library, see #11763.

@arun3688, @AnHeuermann could you please investigate why this happens?

Thanks!

Keeping @deepak19015 and @AndreaBartolini in the loop.

Apr 19 '24 22:04 casella

Further piece of information. If you try this MWE

model TestExponential
  Real x(start = 1, fixed = true);
  Real x_check;
equation
  der(x) = -x;
  x_check = exp(-time);
  assert(abs(x - x_check)<1e-6, "Numerical solution too far from exact one");
  Modelica.Utilities.Streams.print("x = "+String(x,".10f")+", x_check = " + String(x_check, ".10f")+", time = "+String(time));
annotation(
    experiment(StartTime = 0, StopTime = 2, Tolerance = 1e-08, Interval = 0.5));
end TestExponential;

with

from fmpy import *
fmu = './TestExponential.fmu'
dump(fmu)
result = simulate_fmu(fmu,fmi_type='ModelExchange',output_interval=0.5)

you get

LOG_STDOUT        | info    | x = 1.0000000000, x_check = 1.0000000000, time = 0
LOG_STDOUT        | info    | x = 1.0000000000, x_check = 1.0000000000, time = 0
LOG_STDOUT        | info    | x = 1.0000000000, x_check = 1.0000000000, time = 0
LOG_STDOUT        | info    | x = 0.6065306106, x_check = 0.6065306597, time = 0.5
LOG_STDOUT        | info    | x = 0.6065306106, x_check = 0.6065306597, time = 0.5
LOG_STDOUT        | info    | x = 0.3678794141, x_check = 0.3678794412, time = 1
LOG_STDOUT        | info    | x = 0.3678794141, x_check = 0.3678794412, time = 1
LOG_STDOUT        | info    | x = 0.2231301458, x_check = 0.2231301601, time = 1.5
LOG_STDOUT        | info    | x = 0.2231301458, x_check = 0.2231301601, time = 1.5
LOG_STDOUT        | info    | x = 0.1353352758, x_check = 0.1353352832, time = 2
LOG_STDOUT        | info    | x = 0.1353352758, x_check = 0.1353352832, time = 2

In this case, the print statement is still issued twice per time step (which is unnecessary), but at least the result is OK both times. I'm not really sure what is the difference between this MWE and the previous one, but this could also help.

Apr 19 '24 23:04 casella