Open mwetter opened 3 years ago
@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.
@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 bfc5443cb
Todo:
AirHeater_T
) and/or fail to simulate (Resources/Scripts/Dymola/Fluid/HeatPumps/Calibration/ScrollWaterToWater.mos
)@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.
@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.
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.
That's a fair point. Chances are that the solver would still pick T
as a state, though.
@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 https://github.com/lbl-srg/modelica-buildings/issues/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
@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?
It was T(nominal=100)
. I will try a variant with dT=T-reference_T
, and then we can set it to 10.
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 096aa2adc2f5fd95c63785785043c576778a2e06 and Dymola 2020x. The trajectory of the above command is shown below:
Todo:
The latest development branches are
issue1412_stateSelect_Tp
andissue1412_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.
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:
In the Buildings library result comparison across tools, we remove the inverse function annotation in
basicFlowFunction_m_flow
and inbasicFlowFunction_dp
. This leads toIt 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.