Change selection of state variables

[mwetter]

The model IBPSA.Fluid.Interfaces.Examples.Humidifier_u shows oscillations if simulated with CVODE:

If changed to Radau, these oscillations do not occur.

Moreover, the statistics is for the original formulation:

Sizes of linear systems of equations: {3, 3, 3, 3}
Sizes after manipulation of the linear systems: {0, 0, 0, 0}
Sizes of nonlinear systems of equations: {3, 3, 1, 4, 3, 3, 4, 1, 4, 3}
Sizes after manipulation of the nonlinear systems: {1, 1, 0, 1, 1, 1, 1, 0, 1, 1}
Number of numerical Jacobians: 0

Initialization problem
Sizes of nonlinear systems of equations: {3, 1, 3, 1, 4, 1, 3, 1, 3, 1, 4, 1, 4, 1, 3, 1}
Sizes after manipulation of the nonlinear systems: {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}
Number of numerical Jacobians: 0

In the Buildings library result comparison across tools, we remove the inverse function annotation in basicFlowFunction_m_flow and in basicFlowFunction_dp. This leads to

Sizes of linear systems of equations: {3, 3, 3, 3}
Sizes after manipulation of the linear systems: {0, 0, 0, 0}
Sizes of nonlinear systems of equations: {3, 3, 1, 4, 3, 3, 4, 1, 4, 3}
Sizes after manipulation of the nonlinear systems: {1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
Number of numerical Jacobians: 0

Initialization problem
Sizes of nonlinear systems of equations: {3, 3, 1, 4, 3, 3, 4, 1, 4, 3}
Sizes after manipulation of the nonlinear systems: {1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
Number of numerical Jacobians: 0

It turns out that this model then does not simulate with Cvode 1E-6 because the oscillations trigger the assertions in the media, e.g., the model is unstable. Best is to switch for this model to the Radau solver.

[Mathadon]

@mwetter this is a workaround while the root cause of the problem is inadequate tolerance settings. The state variable mix2.mXi[1] has a nominal value of 1e-8 for the given model parameters, which the time integrator solves to too low accuracy, causing the oscillations. Adding the nominal attribute to this variable seems to solve the problem too:

  Modelica.SIunits.Mass[Medium.nXi] mXi(
    nominal=fluidVolume*1000*X_start[1:Medium.nXi],
    start=fluidVolume*rho_start*X_start[1:Medium.nXi])

Note that the start attribute was already there, but it is incorrect since the density of the main fluid is used instead of the trace substance. Therefore I hardcoded a density of 1000 kg/m3 for water. A better solution could be to grab it from the Medium properties but it does not seem to exist and it certainly is not defined for all media. 1000kg/s seems to be a good default though.

I propose to revise the PR #1413 accordingly.

[mwetter]

@Mathadon : A problem is that fluidVolume is not know during compilation (non-literal value) and hence the nominal attribute will be ignored. I tested another implementation that instead uses medium.Xi as the state so we don't need to care about fluidVolume*density (density should be 1.2, not 1000). Then I still get the oscillations with CVode. This is commit bfc5443

[mwetter]

Todo:

• check models that have larger systems of nonlinear equations (e.g., AirHeater_T) and/or fail to simulate (Resources/Scripts/Dymola/Fluid/HeatPumps/Calibration/ScrollWaterToWater.mos)
• Compare computing time in Buildings library for large models

[mwetter]

@Mathadon : I did some further modifications (all pushed to git) to make u a state also for water, and to avoid state select if the ConservationEquation is used for a steady-state model.
Some models have now more non-linear equations. I am not sure if this is only the case when volumes are connected to volumes in which case an index reduction is needed. I have not looked into whether these can be avoided.

Also, if we use mSenFac > 1, then T becomes a state rather than u. In the past, we had always U as a state, see for example Buildings.Fluid.MixingVolumes.Validation.MixingVolumeMFactor.

In short, it looks like U as a state gives better equations, but the nominal value may not be a literal value. I wonder if we should end up looking into introducing _u = U/(V*rho_start) in ConservationEquation, and then have _u as a state. Then everything is known.

[Mathadon]

@mwetter
As far as I can tell, index reduction on mixing volumes should cause linear algebraic loops since there is a linear relationship between mass and pressure. So the non-linear equations probably are related to something else. What models show these different statistics?

For mSenFac>1, it makes sense that T is a state since our equations allow computing u and U from T:

U = m*medium.u + CSen*(medium.T-Medium.reference_T);
//baseproperties:
    dT = T - reference_T;
    h = dT*dryair.cp * X_air +
       (dT * steam.cp + h_fg) * X_steam;
    u = h-pStp/dStp;

I think that we'll have to do some more equation inverting (in the equations above) if we want u to be selected as a state.

[mwetter]

I think Buildings.Fluid.MixingVolumes.Validation.MixingVolumeMFactor is a good example.

What do you think about _u = U/(V*rho_start), as this is what the solver is anyway doing with the nominal attribute.

[Mathadon]

That's a fair point. Chances are that the solver would still pick T as a state, though.

[mwetter]

@Mathadon : I run experiments on all models of the Buildings library, using branches IBPSASync_issue1412_humidifer_u_Cvode (ec06bb573182d9952dad3f71348c3b5ac8efd85c) and IBPSASync_issue1412_humidifer_u_Cvode_T_state (c5ab735f774adf0c563e82d123e84692efc06d5d).
The first uses _u as a state, and the second uses T as a state.
These comparisons use only one simulation run for each model, so there is a spread in the results. However, I run some multiple times and the trend stays the same.

Using _u rather than U (the original implementation) leads generally to slightly slower simulation, but I attribute this to having now an actual error control on _u.

For Optimica, the differences are

For Dymola, the differences are

Selecting T rather than _u as a state generally leads to faster simulation with Dymola and OPTIMICA, typically because of fewer Jacobian evaluations. The overhead of computing T rather than directly using _u is not noticeable.

For Optimica, the differences are

For Dymola, the differences are

Therefore, using T rather than _u as a state seems fine.

Another issue I see is that the state is pressure of density. (For Buildings.Fluid.MixingVolumes.Examples.MixingVolumeMoistAir, OCT select pressure while Dymola selects density as a state on branch IBPSASync_issue1412_humidifer_u_Cvode_T_state.)

I see either problematic because they only vary slightly, for example, pressure is 1E5 and variations can be in the order of 1 to 1000. I will try on branch IBPSASync_issue1412_humidifer_u_Cvode_dp_T_state to use T and a relative pressure as the state. See also lbl-srg/modelica-buildings#1671

Edit:
Better use normalized density as a state as density appears in equations such as the one below from DualFanDualDuct:

fanRet.vol.dynBal.medium._der_d := fanRet.vol.dynBal.mb_flow/ fanRet.vol.V

Branch will be IBPSASync_issue1412_humidifer_u_Cvode_d_T_state

[Mathadon]

@mwetter what is the nominal value for T in these simulations when T is a state? Nominal temperature differences are in the order of 1-2K for zone air while the nominal value of T is about 290K. This could cause the solver to solve for lower accuracy when T is a state. I guess that these benchmarks don't take solver accuracy into account?

[mwetter]

It was T(nominal=100). I will try a variant with dT=T-reference_T, and then we can set it to 10.

[mwetter]

Note that with the old state selection,

simulateModel("IBPSA.Airflow.Multizone.Validation.OpenDoorBuoyancyPressureDynamic", stopTime=14400, method="Cvode", tolerance=1e-06, resultFile="OpenDoorBuoyancyPressureDynamic");

simulates very slow. With Radau 1E-6, and with CVode 1E-8, it simulates quickly. This is with 096aa2a and Dymola 2020x. The trajectory of the above command is shown below:

[mwetter]

Todo:

• Rerun experiments
• Check which variables appear in nonlinear equation.

[mwetter]

The latest development branches are

• issue1412_stateSelect_Tp and
• issue1412_stateSelect_Tp_offset. This is the same as above, but for air, T rather than T_degC is used as a state, and dp has 1000 Pa added.

The respective branches in Buildings have the same name, but IBPSASync_ added.