Convergence issues on nonlinear equations not involving der()

[Francesco Casella]

Please consider the ThermoPower.Test.WaterComponents.TestValves. It takes 11 seconds to simulate in OMC and 0.12 seconds in Dymola, with +cseCall active.

Turning on LOG_NLS_V reveals the following behaviour:

• at time 3.0000 the NLS requires 137 iterations to converge, which probably points to badly chosen scaling factors or tolerances
• the main solver takes a step to 3.008
• the NLS now tries for 147 iterations, then it keeps trying doing other things for a long time, then eventually gives up
• the main solver steps back to 3.00000004, NLS converges quickly
• the main solver steps to 3.008, this time the convergence is very fast
• the simulation proceeds normally

Bottom line 1: at time 3.000 there are probably too many iterations.

Bottom line 2: the first attempt at time 3.008 should have probably been aborted earlier, retrying a shorter step.

[Martin Sjölund]

Milestone changed to 1.9.3 since 1.9.2 was released.

[Francesco Casella]

I have managed to perform some further analysis. The simulated transient is such that the pressure V5.fluidState.p drops very quickly from about 7e5 to just 1e5 between time = 2.96 and time = 3.00. This is an excerpt of the info taken from LOG_NLS_V:

• at time = 2.984, the noninear solver takes 4 iterations to converge to p = 410239.35
• at time = 2.992, the nonlinear solver takes 5 iterations to converge to p = 209937.07
• at time = 3.000, the solver extrapolates p = 9634.7855 (which is an excessively low value!) due to the very steep change of this variable over the previous time step. The actual solution at this point in time is p = 100000.00.
• Iteration 1: p = 190020.75
• Iteration 2: p = 10341.522
• Iteration 3: p = 189311.11
• Iteration 4: p = 11053.82
• Iteration 5: p = 188595.88
• and so on and so forth, until convergence at Iteration 136 (!)

It is true this problem is very severely nonlinear, but on the other hand I would expect that a proper damping strategy should rapidly kill these oscillations and lead to a quick convergence. I understand this high number of iterations is due to largely insufficient damping.

Can you see if the damping strategy can be improved to deal with this case more effectively?

[Francesco Casella]

Turning off removeSimpleEquations has a beneficial effect on the convergence of this system: although there are now a few more tearing variables cseNN with equations $cseNN = max(flowVarMM,eps), the solver only undergoes 6 iterations. The first two are the same as before, but then the solver rapidly converges on the target.

• Iteration 1: p = 190020.74
• Iteration 2: p = 10341.531
• Iteration 3: p = 100257.38
• Iteration 4: p = 99999.468
• Iteration 5: p = 100000
• Iteration 6: p = 100000

I attach the output of dumpSimCode in the two cases for your reference.

This is what happens at time = 3.000. The valves upstream V5 are completely closed, so that the flow through V5 (which is open) becomes zero, and as a consequence its dp becomes zero, so that its inlet pressure becomes equal to the outlet pressure, that is 100000. Since the valve equation is explicit in the flow rate, this requires inverting the m_flow = f(dp) function, which is close to a square root and highly nonlinear.

In the case with removeSimpleEquations (which has trouble converging), these are the crucial equations:

79: V5.dp=V5.fluidState.p - SinkP1.p0
83: V3.w= <expression returning zero because V3 is closed>
89: V2.w= <expression returning zero because V2 is closed>
90: V5.w=V2.w + V3.w // that is, zero
102: V5.rho=Modelica.Media.Water.IF97_Utilities.rho_props_ph(V5.fluidState.p, DIVISION($cse1 * V2.outlet.h_outflow + $cse2 * V3.outlet.h_outflow, $cse1 + $cse2), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(V5.fluidState.p, DIVISION($cse1 * V2.outlet.h_outflow + $cse2 * V3.outlet.h_outflow, $cse1 + $cse2), 0, 0))[Real ]
106: homotopy(ThermoPower.Water.ValveLiq$V5.FlowChar(CloseLoad.y) * V5.Av * sqrt(V5.rho) * ThermoPower.Water.ValveLiq$V5.sqrtR(V5.dp, 4000.0), CloseLoad.y * V5.wnom * DIVISION(V5.fluidState.p - SinkP1.p0, V5.thetanom * V5.dpnom)) - V5.w (RESIDUAL)

Considering lambda = 1 and the values of previously computed variables, this boils down to

V5.dp=V5.fluidState.p - SinkP1.p0
V5.w=0
V5.rho=IF97_density(V5.fluidState.p, some_mixing_enthalpy)
V5.Av * sqrt(V5.rho) * ThermoPower.Water.ValveLiq$V5.sqrtR(V5.dp, 4000.0) - V5.w (RESIDUAL)

In the case without removeSimpleEquation, we get

236: V5.dp=V5.inlet.p - V5.outlet.p
237: $cse36=ThermoPower.Water.ValveLiq$V5.sqrtR(V5.dp, 4000.0)
245: $cse39=Modelica.Media.Water.IF97_Utilities.rho_props_ph(V5.inlet.p, DIVISION($cse26 * V2.outlet.h_outflow + $cse17 * V3.outlet.h_outflow, $cse26 + $cse17), Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(V5.inlet.p, DIVISION($cse26 * V2.outlet.h_outflow + $cse17 * V3.outlet.h_outflow, $cse26 + $cse17), 0, 0))
246: V5.fluidState.d=$cse39
247: V5.rho=V5.fluidState.d
248: $cse35=sqrt(V5.rho)
249: V5.w=$cse34 * V5.Av * $cse35 * $cse36
250: V5.inlet.m_flow=V5.w
various equations leading to V3.outlet.m_flow = 0, V2.outlet.m_flow = 0
289: V3.outlet.m_flow + V5.inlet.m_flow + V2.outlet.m_flow (RESIDUAL)

As far as I understand, apart from involving a bit more variables because of not removing aliases, the process that leads to the computation of the inlet pressure of V5 is the same: the pressure is the tearing variable, from the pressure you compute the density, then the flow, then the flow is brought to zero by Newton iterations. In this case, the variability of enthalpies plays a marginal role.

The residual equation mentioned above is very strongly nonlinear for two reasons. One is the dependency of the flow on the (regularized) square root of dp. As the square root is a concave function, this typically leads to large overshoots and oscillations of Newton-type methods which are not properly damped.

The other reason is that the flow depends on the square root of the density. As long as the pressure is greater than the saturation pressure (which is about 3000 in this case), sqrt(rho) is virtually constant. It seems that this is the case here, so the only nonlinearity that the Newton solver has to fight against is that of the sqrtR() function.

Summing up, with default options the convergence is very slow and highly oscillatory. Without removeSimpleEquations, virtually the same system is obtained. The first two iterations give the same results, but for some reason the third iteration gets right on target, instead of continuing the oscillations.

Please investigate the reason of this difference thoroughly.

[Francesco Casella]

Dump of simcode

[Willi Braun]

Replying to casella:

Can you see if the damping strategy can be improved to deal with this case more effectively?

You are right currently the damping strategy fails for this example completely, after playing a bit with some constants it's possible to bring this 136 iteration to 6. Now I need to cross check that with other examples.

[Willi Braun]

Fixed damping strategy in fff5d10/OMCompiler.

[Francesco Casella]

I still don't see any improvement in the testsuite: in both

​https://test.openmodelica.org/libraries/ThermoPower/BuildModelRecursive.html

and

​https://test.openmodelica.org/libraries/ThermoPower_Experimental/BuildModelRecursive.html

(which do not have removeSimpleEquations), the simulation time for the test is still reported at 7 seconds for simulation, as in the past. I am currently unable to check thoroughly the LOG_NLS_V output with a local installation, but I suspect the problem is still there, otherwise the simulation time should have been reduced significantly.

[Martin Sjölund]

Moved to new milestone 1.9.4

[Martin Sjölund]

Milestone pushed to 1.9.5

[Martin Sjölund]

Milestone renamed

[Martin Sjölund]

Ticket retargeted after milestone closed

[Martin Sjölund]

Milestone moved to 1.12.0 due to 1.11.0 already being released.

[Francesco Casella]

Moved to milestone 1.13.0 as milestone 1.12.0 was closed

[Francesco Casella]

Milestone moved to 1.13.0 due to 1.12.0 already being released.

[Francesco Casella]

Rescheduled to 1.14.0 after 1.13.0 releasee

[Francesco Casella]

Releasing 1.14.0 which is stable and has many improvements w.r.t. 1.13.2. This issue is rescheduled to 1.16.0

[Francesco Casella]

Retargeted to 1.17.0 after 1.16.0 release

[Francesco Casella]

Retargeted to 1.18.0 because of 1.17.0 timed release.

[Francesco Casella]

Ticket retargeted after milestone closed