id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
5005	Issue with DASSL/IDA and implicit systems with non-smooth functions	Francesco Casella	Willi Braun	"Please check the results of the verification of the [https://libraries.openmodelica.org/branches/master/Modelica_trunk/files/Modelica_trunk_Modelica.Electrical.Analog.Examples.AmplifierWithOpAmpDetailed.diff.html Modelica.Electrical.Analog.Examples.AmplifierWithOpAmpDetailed]. Starting from time = 0.0008946, the results start diverging, see, e.g, [https://libraries.openmodelica.org/branches/master/Modelica_trunk/files/Modelica_trunk_Modelica.Electrical.Analog.Examples.AmplifierWithOpAmpDetailed.diff.opAmp.v_in.html opAmp.v_in].

I made some analysis of what happens there and I report the conclusions here. If you use the declarative debugger, it turns out that all system equations are analytic execept for the following three ones:
{{{
(assign) der(opAmp.v_source) := smooth(0, if noEvent(der(opAmp.x) > opAmp.sr_p_val) then opAmp.sr_p_val else if noEvent(der(opAmp.x) < opAmp.sr_m_val) then opAmp.sr_m_val else der(opAmp.x))
(assign) opAmp.q_sum_help := Modelica.Electrical.Analog.Basic.OpAmpDetailed.FCNq_sum_limit(opAmp.q_sum, opAmp.q_fp1, constantVoltage.V, constantVoltage1.V, opAmp.vcp_abs, opAmp.vcm_abs)
(torn) opAmp.i_out := Modelica.Electrical.Analog.Basic.OpAmpDetailed.FCNiout_limit(opAmp.v_source, resistor2.v, opAmp.Rout, opAmp.Imaxsi_val, opAmp.Imaxso_val)
}}}
The first two are assignments in the causalized system, while the third shows up in a nonlinear torn implicit system of equations, using {{{resistor2.v}}} as a tearing variable.

As far as I understand, in all three cases the right-hand sides are continuous functions, but not smooth ones, as hard limiters are involved. Since noEvent() and non-inlined functions are involved, this strains the solvers a bit. In fact, I'm not sure why the decision was made to avoid events in this case, but I understand the idea is that as long as there are no hard discontinuities, solvers should manage without events. BTW, Dymola (which was used to generate the reference results) does manage without any issue.

I tried to replicate the results reported on Hudson on my local OMEdit setup. If I run the example out of the box, I get the correct result up to time = 0.0013944, then it starts to deviate from the reference (and correct) one. This is possibly because some of the settings on Hudson are different, I'm not sure.

I then tried to change the ODE  solver to IDA, and/or the nonlinear solver to Kinsol, and/or to switch off Equidistant time grid. Each combination of these factors causes the simulation to veer off course at different points, but none of them provided the correct results. Also {{{--daeMode=new}}} leads to the same kind of problems.

The only way I could get the correct result was to set {{{--maxSizeNonlinearTearing=0}}}. DASSL provided the correct result right away, while IDA still got some glitches and required to set {{{-noEquidistantTimeGrid}}} to get the correct result. IDA in {{{--daeMode=new}}} still failed to produce the correct results.

I guess this means that the interaction between the ODE/DAE solver and the nonlinear equations with non-smooth terms is to blame here. BTW, Dymola selects a different tearing variable, {{{resistor2.i}}}, which may explain why it works fine.

So, my conclusion is that the way we currently deal with these systems, that contain smooth(0, ..) expressions and/or continuous but non-smooth functions, is numerically fragile, for reasons that need to be further investigated.  

I believe this is a serious issue that probably affects many other models in the MSL and in other libraries, so it should be fixed with high priority, and most definitely addressed for version 2.0.0"	defect	assigned	blocker	2.0.0	Run-time				Bernhard Bachmann