Context Navigation

← Previous Ticket
Next Ticket →

#5770 assigned defect

The backend incorrectly generates a large set of nonlinear equations

Reported by:	Francesco Casella	Owned by:	Karim Adbdelhak
Priority:	high	Milestone:
Component:	Backend	Version:
Keywords:		Cc:	Andreas Heuermann, giovanni.mangola@…, matthis.thorade@…

Description

@Karim, I managed to come up with a test case for the slow behaviour in the power plant simulation that I mentioned during today's webmeeting. The test case has a reduced size to be more manageable.

Please run the attached Save Total model. The backend identifies a nonlinear system with 18 unknowns and 7 iteration variables:

der(fWPH_11_1.flowSTEAM.htilde[1])
der(fWPH_11_1.flowSTEAM.htilde[2])
der(fWPH_11_1.nusseltCondenser.hl)
der(fWPH_11_1.nusseltCondenser.hv)
der(fWPH_11_1.nusseltCondenser.p)
der(fWPH_11_1.nusseltCondenser.Vv)
fWPH_11_1.nusseltCondenser.ws_in

Dymola identifies a linear system instead, which has the 6 following unknowns after tearing.

der(fWPH_11_1.flowSTEAM.h[2])
der(fWPH_11_1.flowSTEAM.h[3])
der(fWPH_11_1.nusseltCondenser.hl)
der(fWPH_11_1.nusseltCondenser.hv)
der(fWPH_11_1.nusseltCondenser.p)
der(fWPH_11_1.nusseltCondenser.Vl)

The first two and the last ones aliases of the first two and of the sixth found by OMC.

The simulation of this model is actually a steady-state, so there are only 19 variable step-size steps, and 500 communication intervals. Yet, the nonlinear system is reported by the debugger to be executed a whopping 48209 times. In the case of larger systems, this can be really become the bottleneck.

I scanned the equations of the nonlinear system in the OMC debugger, but I couldn't understand why OMC considers them nonlinear. Could you please check?

Attachments (2)

TestTotal.mo (454.9 KB ) - added by Francesco Casella 5 years ago.
test.mos (136 bytes ) - added by Francesco Casella 5 years ago.

Download all attachments as: .zip

Change History (31)

by Francesco Casella, 5 years ago

Attachment:	TestTotal.mo added

comment:1 by Francesco Casella, 5 years ago

Cc:	Andreas Heuermann added; Karim Adbdelhak removed
Owner:	changed from Lennart Ochel to Karim Adbdelhak
Status:	new → assigned

comment:2 by Francesco Casella, 5 years ago

Cc:	giovanni.mangola@… added

comment:3 by Andreas Heuermann, 5 years ago

What are the needed flags to run this example? I'm trying with this runTest.mos

loadFile("TestTotal.mo"); getErrorString();
setComandLineOptions("-d=newInst"); getErrrorString();
simulate(Test); getErrorString();

but get only an error:

$ omc runTest.mos
true
""
record SimulationResult
    resultFile = "",
    simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'Test', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
    messages = "Failed to build model: Test",
    timeFrontend = 10.083323602,
    timeBackend = 0.0,
    timeSimCode = 0.0,
    timeTemplates = 0.0,
    timeCompile = 0.0,
    timeSimulation = 0.0,
    timeTotal = 10.08335457
end SimulationResult;
"[/home/aheuermann1/workspace/Testitesttest/ticket/ticket5770/TestTotal.mo:348:354-348:381:writable] Notification: From here:
[/home/aheuermann1/workspace/Testitesttest/ticket/ticket5770/TestTotal.mo:1327:7-1327:91:writable] Error: Array dimension mismatch, expression {fWPH_11_1.Tstart_wall_des[1]} has type Real(start = 500.0, nominal = 500.0, min = 273.15, max = 2273.15, quantity = "ThermodynamicTemperature", unit = "K", displayUnit = "degC")[1], expected array dimensions [2].
Error: Error occurred while flattening model Test
"

comment:4 by Francesco Casella, 5 years ago

I guess the reasons are the typos in your script of comment:3.

I attach the test script, sorry for not doing it earlier :)

by Francesco Casella, 5 years ago

Attachment:	test.mos added

follow-up: 19 comment:5 by Karim Adbdelhak, 5 years ago

I updated the -d=dumpLoops option for better visual partitioning of components and added -d=dumpLoopsVerbose to dump solvability information about the components in PR664. Here is a sample output for the component in question:

################################################################################
 dumpLoops: SORTED COMPONENT 
################################################################################

torn nonlinear Equationsystem:

internal vars
1: $DER.fWPH_11_1.nusseltCondenser.El:DUMMY_DER(fixed = false )  "Condensate liquid energy" type: Real
2: $DER.fWPH_11_1.nusseltCondenser.Ev:DUMMY_DER(fixed = false )  "Vapour energy" type: Real
3: $DER.fWPH_11_1.nusseltCondenser.Mv:DUMMY_DER(fixed = false )  "Vapour mass" type: Real
4: $DER.fWPH_11_1.nusseltCondenser.Ml:DUMMY_DER(fixed = false )  "Condensate liquid mass" type: Real
5: $DER.fWPH_11_1.nusseltCondenser.rho_l:DUMMY_DER(fixed = false )  "Liquid density" type: Real
6: $DER.fWPH_11_1.nusseltCondenser.rho_v:DUMMY_DER(fixed = false )  "Steam density" type: Real
7: $DER.fWPH_11_1.flowSTEAM.p:DUMMY_DER(fixed = false )  "Fluid pressure for property calculations" type: Real
8: fWPH_11_1.flowSTEAM.dMdt[2]:VARIABLE(unit = "kg/s" protected = true )  "Derivative of fluid mass in each volume" type: Real [2]
9: fWPH_11_1.flowSTEAM.dMdt[1]:VARIABLE(unit = "kg/s" protected = true )  "Derivative of fluid mass in each volume" type: Real [2]
10: fWPH_11_1.flowSTEAM.wbar[2]:VARIABLE(min = -100000.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom / /*Real*/(fWPH_11_1.flowSTEAM.Nt) unit = "kg/s" )  "Average mass flow rates (single tube)" type: Real [2]
11: fWPH_11_1.flowSTEAM.wbar[1]:VARIABLE(min = -100000.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom / /*Real*/(fWPH_11_1.flowSTEAM.Nt) unit = "kg/s" )  "Average mass flow rates (single tube)" type: Real [2]

residual vars
1: fWPH_11_1.nusseltCondenser.p:STATE(1)(min = 611.657 max = 100000000.0 start = fWPH_11_1.flowSTEAM.pstart unit = "Pa" nominal = 1000000.0 stateSelect=StateSelect.prefer )  "Medium pressure" type: Real
2: fWPH_11_1.nusseltCondenser.ws_in:VARIABLE(min = 0.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom unit = "kg/s" )  "Steam mass flow rate" type: Real
3: fWPH_11_1.nusseltCondenser.hl:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.nusseltCondenser.hl_start unit = "J/kg" nominal = 500000.0 stateSelect=StateSelect.prefer )  "Specific liquid enthalpy" type: Real
4: $DER.fWPH_11_1.nusseltCondenser.Vv:DUMMY_DER(fixed = false )  "Volume occupied by vapour" type: Real
5: fWPH_11_1.nusseltCondenser.hv:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.nusseltCondenser.hv_start unit = "J/kg" nominal = 500000.0 stateSelect=StateSelect.prefer )  "Specific vapor enthalpy" type: Real
6: fWPH_11_1.flowSTEAM.htilde[2]:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.flowSTEAM.hstart[3] unit = "J/kg" nominal = 500000.0 )  "Enthalpy state variables" type: Real [2]
7: fWPH_11_1.flowSTEAM.htilde[1]:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.flowSTEAM.hstart[2] unit = "J/kg" nominal = 500000.0 )  "Enthalpy state variables" type: Real [2]

internal equation
1/1 (1): $DER.fWPH_11_1.nusseltCondenser.El = fWPH_11_1.nusseltCondenser.ws_in * (1.0 - fWPH_11_1.nusseltCondenser.x_s) * $cse51 + DRAINin.w0 * (1.0 - fWPH_11_1.nusseltCondenser.x_d) * $cse52 + fWPH_11_1.nusseltCondenser.Q_flux + (-fWPH_11_1.nusseltCondenser.wd_out) * fWPH_11_1.nusseltCondenser.hl - fWPH_11_1.nusseltCondenser.w_ev * fWPH_11_1.nusseltCondenser.hv_sat   [dynamic |0|0|0|0|]
2/2 (1): $DER.fWPH_11_1.nusseltCondenser.Ev = fWPH_11_1.nusseltCondenser.ws_in * fWPH_11_1.nusseltCondenser.x_s * $cse53 + DRAINin.w0 * fWPH_11_1.nusseltCondenser.x_d * $cse54 + fWPH_11_1.nusseltCondenser.w_ev * fWPH_11_1.nusseltCondenser.hv_sat + (-fWPH_11_1.nusseltCondenser.Q_cond) - fWPH_11_1.nusseltCondenser.Q_flux   [dynamic |0|0|0|0|]
3/3 (1): $DER.fWPH_11_1.nusseltCondenser.Mv = fWPH_11_1.nusseltCondenser.ws_in * fWPH_11_1.nusseltCondenser.x_s + DRAINin.w0 * fWPH_11_1.nusseltCondenser.x_d + fWPH_11_1.nusseltCondenser.w_ev - fWPH_11_1.nusseltCondenser.w_cond   [dynamic |0|0|0|0|]
4/4 (1): $DER.fWPH_11_1.nusseltCondenser.Ml = fWPH_11_1.nusseltCondenser.ws_in * (1.0 - fWPH_11_1.nusseltCondenser.x_s) + DRAINin.w0 * (1.0 - fWPH_11_1.nusseltCondenser.x_d) + fWPH_11_1.nusseltCondenser.w_cond + (-fWPH_11_1.nusseltCondenser.wd_out) - fWPH_11_1.nusseltCondenser.w_ev   [dynamic |0|0|0|0|]
5/5 (1): $DER.fWPH_11_1.nusseltCondenser.rho_l = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hl, $cse47, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hl))   [unknown |0|0|0|0|]
6/6 (1): $DER.fWPH_11_1.nusseltCondenser.rho_v = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hv, $cse46, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hv))   [unknown |0|0|0|0|]
7/7 (1): $DER.fWPH_11_1.flowSTEAM.p = der(fWPH_11_1.nusseltCondenser.p)   [dynamic |0|0|0|0|]
8/8 (1): fWPH_11_1.flowSTEAM.dMdt[2] = 4.2 * fWPH_11_1.flowSTEAM.A * (fWPH_11_1.flowSTEAM.drbdh1[2] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.drbdh2[2] * der(fWPH_11_1.flowSTEAM.htilde[2]) + fWPH_11_1.flowSTEAM.drbdp[2] * $DER.fWPH_11_1.flowSTEAM.p)   [dynamic |0|0|0|0|]
9/9 (1): fWPH_11_1.flowSTEAM.dMdt[1] = 4.2 * fWPH_11_1.flowSTEAM.A * (fWPH_11_1.flowSTEAM.drbdh2[1] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.drbdp[1] * $DER.fWPH_11_1.flowSTEAM.p)   [dynamic |0|0|0|0|]
10/10 (1): fWPH_11_1.flowSTEAM.wbar[2] = fWPH_11_1.flowSTEAM.heatTransfer.w[3] + (-0.5) * fWPH_11_1.flowSTEAM.dMdt[2] - fWPH_11_1.flowSTEAM.dMdt[1]   [dynamic |0|0|0|0|]
11/11 (1): fWPH_11_1.flowSTEAM.wbar[1] = fWPH_11_1.flowSTEAM.heatTransfer.w[3] + (-0.5) * fWPH_11_1.flowSTEAM.dMdt[1]   [dynamic |0|0|0|0|]

residual equations
1/1 (1): 4.2 * fWPH_11_1.flowSTEAM.A * fWPH_11_1.flowSTEAM.rho[2] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.wbar[1] * (fWPH_11_1.flowSTEAM.htilde[1] - STEAMin.h) + 4.2 * (-fWPH_11_1.flowSTEAM.A) * $DER.fWPH_11_1.flowSTEAM.p = fWPH_11_1.flowSTEAM.Q_single[1]   [dynamic |0|0|0|0|]
2/2 (1): fWPH_11_1.flowSTEAM.dMdt[1] + fWPH_11_1.flowSTEAM.dMdt[2] = (fWPH_11_1.STEAMin.m_flow - fWPH_11_1.nusseltCondenser.ws_in) / /*Real*/(fWPH_11_1.flowSTEAM.Nt)   [dynamic |0|0|0|0|]
3/3 (1): 4.2 * fWPH_11_1.flowSTEAM.A * fWPH_11_1.flowSTEAM.rhobar[2] * der(fWPH_11_1.flowSTEAM.htilde[2]) + fWPH_11_1.flowSTEAM.wbar[2] * (fWPH_11_1.flowSTEAM.htilde[2] - fWPH_11_1.flowSTEAM.htilde[1]) + 4.2 * (-fWPH_11_1.flowSTEAM.A) * $DER.fWPH_11_1.flowSTEAM.p = fWPH_11_1.flowSTEAM.Q_single[2]   [dynamic |0|0|0|0|]
4/4 (1): $DER.fWPH_11_1.nusseltCondenser.Ev = fWPH_11_1.nusseltCondenser.Mv * (der(fWPH_11_1.nusseltCondenser.hv) + (fWPH_11_1.nusseltCondenser.p * $DER.fWPH_11_1.nusseltCondenser.rho_v - der(fWPH_11_1.nusseltCondenser.p) * fWPH_11_1.nusseltCondenser.rho_v) / fWPH_11_1.nusseltCondenser.rho_v ^ 2.0) + $DER.fWPH_11_1.nusseltCondenser.Mv * (fWPH_11_1.nusseltCondenser.hv - fWPH_11_1.nusseltCondenser.p / fWPH_11_1.nusseltCondenser.rho_v)   [dynamic |0|0|0|0|]
5/5 (1): $DER.fWPH_11_1.nusseltCondenser.Mv = fWPH_11_1.nusseltCondenser.Vv * $DER.fWPH_11_1.nusseltCondenser.rho_v + $DER.fWPH_11_1.nusseltCondenser.Vv * fWPH_11_1.nusseltCondenser.rho_v   [dynamic |0|0|0|0|]
6/6 (1): $DER.fWPH_11_1.nusseltCondenser.Ml = fWPH_11_1.nusseltCondenser.Vl * $DER.fWPH_11_1.nusseltCondenser.rho_l - $DER.fWPH_11_1.nusseltCondenser.Vv * fWPH_11_1.nusseltCondenser.rho_l   [dynamic |0|0|0|0|]
7/7 (1): $DER.fWPH_11_1.nusseltCondenser.El = fWPH_11_1.nusseltCondenser.Ml * (der(fWPH_11_1.nusseltCondenser.hl) + (fWPH_11_1.nusseltCondenser.p * $DER.fWPH_11_1.nusseltCondenser.rho_l - der(fWPH_11_1.nusseltCondenser.p) * fWPH_11_1.nusseltCondenser.rho_l) / fWPH_11_1.nusseltCondenser.rho_l ^ 2.0) + $DER.fWPH_11_1.nusseltCondenser.Ml * (fWPH_11_1.nusseltCondenser.hl - fWPH_11_1.nusseltCondenser.p / fWPH_11_1.nusseltCondenser.rho_l)   [dynamic |0|0|0|0|]


################################################################################
 dumpLoopsVerbose: UNSORTED COMPONENT WITH ENHANCED ADJACENCY MATRIX 
################################################################################


component variables (18)
========================================
1: fWPH_11_1.flowSTEAM.htilde[1]:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.flowSTEAM.hstart[2] unit = "J/kg" nominal = 500000.0 )  "Enthalpy state variables" type: Real [2]
2: fWPH_11_1.flowSTEAM.htilde[2]:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.flowSTEAM.hstart[3] unit = "J/kg" nominal = 500000.0 )  "Enthalpy state variables" type: Real [2]
3: fWPH_11_1.nusseltCondenser.hv:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.nusseltCondenser.hv_start unit = "J/kg" nominal = 500000.0 stateSelect=StateSelect.prefer )  "Specific vapor enthalpy" type: Real
4: $DER.fWPH_11_1.nusseltCondenser.Vv:DUMMY_DER(fixed = false )  "Volume occupied by vapour" type: Real
5: fWPH_11_1.nusseltCondenser.hl:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.nusseltCondenser.hl_start unit = "J/kg" nominal = 500000.0 stateSelect=StateSelect.prefer )  "Specific liquid enthalpy" type: Real
6: fWPH_11_1.nusseltCondenser.ws_in:VARIABLE(min = 0.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom unit = "kg/s" )  "Steam mass flow rate" type: Real
7: fWPH_11_1.nusseltCondenser.p:STATE(1)(min = 611.657 max = 100000000.0 start = fWPH_11_1.flowSTEAM.pstart unit = "Pa" nominal = 1000000.0 stateSelect=StateSelect.prefer )  "Medium pressure" type: Real
8: fWPH_11_1.flowSTEAM.wbar[1]:VARIABLE(min = -100000.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom / /*Real*/(fWPH_11_1.flowSTEAM.Nt) unit = "kg/s" )  "Average mass flow rates (single tube)" type: Real [2]
9: fWPH_11_1.flowSTEAM.wbar[2]:VARIABLE(min = -100000.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom / /*Real*/(fWPH_11_1.flowSTEAM.Nt) unit = "kg/s" )  "Average mass flow rates (single tube)" type: Real [2]
10: fWPH_11_1.flowSTEAM.dMdt[1]:VARIABLE(unit = "kg/s" protected = true )  "Derivative of fluid mass in each volume" type: Real [2]
11: fWPH_11_1.flowSTEAM.dMdt[2]:VARIABLE(unit = "kg/s" protected = true )  "Derivative of fluid mass in each volume" type: Real [2]
12: $DER.fWPH_11_1.flowSTEAM.p:DUMMY_DER(fixed = false )  "Fluid pressure for property calculations" type: Real
13: $DER.fWPH_11_1.nusseltCondenser.rho_v:DUMMY_DER(fixed = false )  "Steam density" type: Real
14: $DER.fWPH_11_1.nusseltCondenser.rho_l:DUMMY_DER(fixed = false )  "Liquid density" type: Real
15: $DER.fWPH_11_1.nusseltCondenser.Ml:DUMMY_DER(fixed = false )  "Condensate liquid mass" type: Real
16: $DER.fWPH_11_1.nusseltCondenser.Mv:DUMMY_DER(fixed = false )  "Vapour mass" type: Real
17: $DER.fWPH_11_1.nusseltCondenser.Ev:DUMMY_DER(fixed = false )  "Vapour energy" type: Real
18: $DER.fWPH_11_1.nusseltCondenser.El:DUMMY_DER(fixed = false )  "Condensate liquid energy" type: Real


component equations (18, 18)
========================================
1/1 (1): $DER.fWPH_11_1.nusseltCondenser.El = fWPH_11_1.nusseltCondenser.ws_in * (1.0 - fWPH_11_1.nusseltCondenser.x_s) * $cse51 + DRAINin.w0 * (1.0 - fWPH_11_1.nusseltCondenser.x_d) * $cse52 + fWPH_11_1.nusseltCondenser.Q_flux + (-fWPH_11_1.nusseltCondenser.wd_out) * fWPH_11_1.nusseltCondenser.hl - fWPH_11_1.nusseltCondenser.w_ev * fWPH_11_1.nusseltCondenser.hv_sat   [dynamic |0|0|0|0|]
2/2 (1): $DER.fWPH_11_1.nusseltCondenser.Ev = fWPH_11_1.nusseltCondenser.ws_in * fWPH_11_1.nusseltCondenser.x_s * $cse53 + DRAINin.w0 * fWPH_11_1.nusseltCondenser.x_d * $cse54 + fWPH_11_1.nusseltCondenser.w_ev * fWPH_11_1.nusseltCondenser.hv_sat + (-fWPH_11_1.nusseltCondenser.Q_cond) - fWPH_11_1.nusseltCondenser.Q_flux   [dynamic |0|0|0|0|]
3/3 (1): $DER.fWPH_11_1.nusseltCondenser.Mv = fWPH_11_1.nusseltCondenser.ws_in * fWPH_11_1.nusseltCondenser.x_s + DRAINin.w0 * fWPH_11_1.nusseltCondenser.x_d + fWPH_11_1.nusseltCondenser.w_ev - fWPH_11_1.nusseltCondenser.w_cond   [dynamic |0|0|0|0|]
4/4 (1): $DER.fWPH_11_1.nusseltCondenser.Ml = fWPH_11_1.nusseltCondenser.ws_in * (1.0 - fWPH_11_1.nusseltCondenser.x_s) + DRAINin.w0 * (1.0 - fWPH_11_1.nusseltCondenser.x_d) + fWPH_11_1.nusseltCondenser.w_cond + (-fWPH_11_1.nusseltCondenser.wd_out) - fWPH_11_1.nusseltCondenser.w_ev   [dynamic |0|0|0|0|]
5/5 (1): $DER.fWPH_11_1.nusseltCondenser.rho_l = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hl, $cse47, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hl))   [unknown |0|0|0|0|]
6/6 (1): $DER.fWPH_11_1.nusseltCondenser.rho_v = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hv, $cse46, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hv))   [unknown |0|0|0|0|]
7/7 (1): $DER.fWPH_11_1.flowSTEAM.p = der(fWPH_11_1.nusseltCondenser.p)   [dynamic |0|0|0|0|]
8/8 (1): fWPH_11_1.flowSTEAM.dMdt[2] = 4.2 * fWPH_11_1.flowSTEAM.A * (fWPH_11_1.flowSTEAM.drbdh1[2] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.drbdh2[2] * der(fWPH_11_1.flowSTEAM.htilde[2]) + fWPH_11_1.flowSTEAM.drbdp[2] * $DER.fWPH_11_1.flowSTEAM.p)   [dynamic |0|0|0|0|]
9/9 (1): fWPH_11_1.flowSTEAM.dMdt[1] = 4.2 * fWPH_11_1.flowSTEAM.A * (fWPH_11_1.flowSTEAM.drbdh2[1] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.drbdp[1] * $DER.fWPH_11_1.flowSTEAM.p)   [dynamic |0|0|0|0|]
10/10 (1): fWPH_11_1.flowSTEAM.wbar[2] = fWPH_11_1.flowSTEAM.heatTransfer.w[3] + (-0.5) * fWPH_11_1.flowSTEAM.dMdt[2] - fWPH_11_1.flowSTEAM.dMdt[1]   [dynamic |0|0|0|0|]
11/11 (1): fWPH_11_1.flowSTEAM.wbar[1] = fWPH_11_1.flowSTEAM.heatTransfer.w[3] + (-0.5) * fWPH_11_1.flowSTEAM.dMdt[1]   [dynamic |0|0|0|0|]
12/12 (1): 4.2 * fWPH_11_1.flowSTEAM.A * fWPH_11_1.flowSTEAM.rho[2] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.wbar[1] * (fWPH_11_1.flowSTEAM.htilde[1] - STEAMin.h) + 4.2 * (-fWPH_11_1.flowSTEAM.A) * $DER.fWPH_11_1.flowSTEAM.p = fWPH_11_1.flowSTEAM.Q_single[1]   [dynamic |0|0|0|0|]
13/13 (1): fWPH_11_1.flowSTEAM.dMdt[1] + fWPH_11_1.flowSTEAM.dMdt[2] = (fWPH_11_1.STEAMin.m_flow - fWPH_11_1.nusseltCondenser.ws_in) / /*Real*/(fWPH_11_1.flowSTEAM.Nt)   [dynamic |0|0|0|0|]
14/14 (1): 4.2 * fWPH_11_1.flowSTEAM.A * fWPH_11_1.flowSTEAM.rhobar[2] * der(fWPH_11_1.flowSTEAM.htilde[2]) + fWPH_11_1.flowSTEAM.wbar[2] * (fWPH_11_1.flowSTEAM.htilde[2] - fWPH_11_1.flowSTEAM.htilde[1]) + 4.2 * (-fWPH_11_1.flowSTEAM.A) * $DER.fWPH_11_1.flowSTEAM.p = fWPH_11_1.flowSTEAM.Q_single[2]   [dynamic |0|0|0|0|]
15/15 (1): $DER.fWPH_11_1.nusseltCondenser.Ev = fWPH_11_1.nusseltCondenser.Mv * (der(fWPH_11_1.nusseltCondenser.hv) + (fWPH_11_1.nusseltCondenser.p * $DER.fWPH_11_1.nusseltCondenser.rho_v - der(fWPH_11_1.nusseltCondenser.p) * fWPH_11_1.nusseltCondenser.rho_v) / fWPH_11_1.nusseltCondenser.rho_v ^ 2.0) + $DER.fWPH_11_1.nusseltCondenser.Mv * (fWPH_11_1.nusseltCondenser.hv - fWPH_11_1.nusseltCondenser.p / fWPH_11_1.nusseltCondenser.rho_v)   [dynamic |0|0|0|0|]
16/16 (1): $DER.fWPH_11_1.nusseltCondenser.Mv = fWPH_11_1.nusseltCondenser.Vv * $DER.fWPH_11_1.nusseltCondenser.rho_v + $DER.fWPH_11_1.nusseltCondenser.Vv * fWPH_11_1.nusseltCondenser.rho_v   [dynamic |0|0|0|0|]
17/17 (1): $DER.fWPH_11_1.nusseltCondenser.Ml = fWPH_11_1.nusseltCondenser.Vl * $DER.fWPH_11_1.nusseltCondenser.rho_l - $DER.fWPH_11_1.nusseltCondenser.Vv * fWPH_11_1.nusseltCondenser.rho_l   [dynamic |0|0|0|0|]
18/18 (1): $DER.fWPH_11_1.nusseltCondenser.El = fWPH_11_1.nusseltCondenser.Ml * (der(fWPH_11_1.nusseltCondenser.hl) + (fWPH_11_1.nusseltCondenser.p * $DER.fWPH_11_1.nusseltCondenser.rho_l - der(fWPH_11_1.nusseltCondenser.p) * fWPH_11_1.nusseltCondenser.rho_l) / fWPH_11_1.nusseltCondenser.rho_l ^ 2.0) + $DER.fWPH_11_1.nusseltCondenser.Ml * (fWPH_11_1.nusseltCondenser.hl - fWPH_11_1.nusseltCondenser.p / fWPH_11_1.nusseltCondenser.rho_l)   [dynamic |0|0|0|0|]

Adjacency Matrix Enhanced (row == equation)
====================================
number of rows: 18
1:(18,solved) (6,variable(true)) (-5,variable(true)) 
2:(17,solved) (6,variable(true)) 
3:(16,solved) (6,variable(true)) 
4:(15,solved) (6,variable(true)) 
5:(14,solved) (-7,unsolvable) (-5,unsolvable) (7,unsolvable) (5,unsolvable) 
6:(13,solved) (-7,unsolvable) (-3,unsolvable) (7,unsolvable) (3,unsolvable) 
7:(12,solved) (7,solved) 
8:(11,solved) (1,variable(true)) (2,variable(true)) (12,variable(true)) 
9:(10,solved) (1,variable(true)) (12,variable(true)) 
10:(9,solved) (11,const(true)) (10,constone) 
11:(8,solved) (10,const(true)) 
12:(1,variable(true)) (8,variable(true)) (-1,variable(true)) (12,param(true)) 
13:(10,constone) (11,constone) (6,param(true)) 
14:(2,variable(true)) (9,variable(true)) (-2,variable(true)) (-1,variable(true)) (12,param(true)) 
15:(17,solved) (3,variable(true)) (-7,variable(true)) (13,variable(true)) (7,variable(true)) (16,variable(true)) (-3,variable(true)) 
16:(16,solved) (13,variable(true)) (4,variable(true)) 
17:(15,solved) (14,variable(true)) (4,variable(true)) 
18:(18,solved) (5,variable(true)) (-7,variable(true)) (14,variable(true)) (7,variable(true)) (15,variable(true)) (-5,variable(true)) 


Transpose Adjacency Matrix Enhanced (row == var)
=====================================
number of rows: 18
1:(-14,variable(true)) (-12,variable(true)) (12,variable(true)) (9,variable(true)) (8,variable(true)) 
2:(-14,variable(true)) (14,variable(true)) (8,variable(true)) 
3:(-15,variable(true)) (15,variable(true)) (6,unsolvable) (-6,unsolvable) 
4:(17,variable(true)) (16,variable(true)) 
5:(-18,variable(true)) (18,variable(true)) (5,unsolvable) (-5,unsolvable) (-1,variable(true)) 
6:(13,param(true)) (4,variable(true)) (3,variable(true)) (2,variable(true)) (1,variable(true)) 
7:(18,variable(true)) (-18,variable(true)) (15,variable(true)) (-15,variable(true)) (7,solved) (6,unsolvable) (-6,unsolvable) (5,unsolvable) (-5,unsolvable) 
8:(12,variable(true)) (11,solved) 
9:(14,variable(true)) (10,solved) 
10:(13,constone) (11,const(true)) (10,constone) (9,solved) 
11:(13,constone) (10,const(true)) (8,solved) 
12:(14,param(true)) (12,param(true)) (9,variable(true)) (8,variable(true)) (7,solved) 
13:(16,variable(true)) (15,variable(true)) (6,solved) 
14:(18,variable(true)) (17,variable(true)) (5,solved) 
15:(18,variable(true)) (17,solved) (4,solved) 
16:(16,solved) (15,variable(true)) (3,solved) 
17:(15,solved) (2,solved) 
18:(18,solved) (1,solved)

The first part shows the result of tearing. (splitting into internal/residuals vars and internal/residual equations). The second (verbose) part shows the solvability information. The indices for each var and equation in these dumps is different due to different ordering. Maybe we can find something odd here.

Last edited 5 years ago by Karim Adbdelhak (previous) (diff)

comment:6 by Karim Adbdelhak, 5 years ago

Just as a reference i will post the Backend structure for solvability:

uniontype Solvability
  record SOLVABILITY_SOLVED "Equation is already solved for the variable" end SOLVABILITY_SOLVED;
  record SOLVABILITY_CONSTONE "Coefficient is equal 1 or -1" end SOLVABILITY_CONSTONE;
  record SOLVABILITY_CONST "Coefficient is constant"
    Boolean b "false if the constant is almost zero (<1e-6)";
  end SOLVABILITY_CONST;
  record SOLVABILITY_PARAMETER "Coefficient contains parameters"
    Boolean b "false if the partial derivative is zero";
  end SOLVABILITY_PARAMETER;
  record SOLVABILITY_LINEAR "Coefficient contains variables, is time varying"
    Boolean b "false if the partial derivative is zero";
  end SOLVABILITY_LINEAR;
  record SOLVABILITY_NONLINEAR "The variable occurs non-linear in the equation." end SOLVABILITY_NONLINEAR;
  record SOLVABILITY_UNSOLVABLE "The variable occurs in the equation, but it is not possible to solve
                     the equation for it." end SOLVABILITY_UNSOLVABLE;
  record SOLVABILITY_SOLVABLE "It is possible to solve the equation for the variable, it is not considered
                     how the variable occurs in the equation." end SOLVABILITY_SOLVABLE;
end Solvability;

comment:7 by Karim Adbdelhak, 5 years ago

Analyzing this system by hand yields a better solution, i guess the one dymola also finds. All following indices relate to the second (verbose) dump. One could solve the equation

13/13 (1): fWPH_11_1.flowSTEAM.dMdt[1] + fWPH_11_1.flowSTEAM.dMdt[2] = (fWPH_11_1.STEAMin.m_flow - fWPH_11_1.nusseltCondenser.ws_in) / /*Real*/(fWPH_11_1.flowSTEAM.Nt)   [dynamic |0|0|0|0|]

for fWPH_11_1.nusseltCondenser.ws_in and make that one an internal variable. EQ 13 depends on VAR 10 and 11 (the other ones are not part of the loop or residuals), which are solved in EQ 8 and 9, which both depend only on VAR 12 (the other ones are not part of the loop or residuals), which is solved in EQ 7. And EQ 7 in the end does not depend on any variables other than der(fWPH_11_1.nusseltCondenser.p) which is a residual variable.

Unfortunately Cellier Tearing greedily selects fWPH_11_1.nusseltCondenser.ws_in as very first variable to be a residual since it appears in so many equations which can directly be causalized by selecting it.

I tried --tearingMethod=omcTearing which is a different approach. It indeed finds the correct set, if one can call it such.

################################################################################
 dumpLoops: SORTED COMPONENT 
################################################################################

torn nonlinear Equationsystem:

internal vars
1: $DER.fWPH_11_1.flowSTEAM.p:DUMMY_DER(fixed = false )  "Fluid pressure for property calculations" type: Real
2: fWPH_11_1.flowSTEAM.dMdt[2]:VARIABLE(unit = "kg/s" protected = true )  "Derivative of fluid mass in each volume" type: Real [2]
3: fWPH_11_1.flowSTEAM.dMdt[1]:VARIABLE(unit = "kg/s" protected = true )  "Derivative of fluid mass in each volume" type: Real [2]
4: fWPH_11_1.flowSTEAM.wbar[2]:VARIABLE(min = -100000.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom / /*Real*/(fWPH_11_1.flowSTEAM.Nt) unit = "kg/s" )  "Average mass flow rates (single tube)" type: Real [2]
5: fWPH_11_1.nusseltCondenser.ws_in:VARIABLE(min = 0.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom unit = "kg/s" )  "Steam mass flow rate" type: Real
6: fWPH_11_1.flowSTEAM.wbar[1]:VARIABLE(min = -100000.0 max = 100000.0 start = fWPH_11_1.flowSTEAM.wnom / /*Real*/(fWPH_11_1.flowSTEAM.Nt) unit = "kg/s" )  "Average mass flow rates (single tube)" type: Real [2]
7: $DER.fWPH_11_1.nusseltCondenser.rho_v:DUMMY_DER(fixed = false )  "Steam density" type: Real
8: $DER.fWPH_11_1.nusseltCondenser.Mv:DUMMY_DER(fixed = false )  "Vapour mass" type: Real
9: $DER.fWPH_11_1.nusseltCondenser.Ev:DUMMY_DER(fixed = false )  "Vapour energy" type: Real
10: $DER.fWPH_11_1.nusseltCondenser.rho_l:DUMMY_DER(fixed = false )  "Liquid density" type: Real
11: $DER.fWPH_11_1.nusseltCondenser.Ml:DUMMY_DER(fixed = false )  "Condensate liquid mass" type: Real
12: $DER.fWPH_11_1.nusseltCondenser.El:DUMMY_DER(fixed = false )  "Condensate liquid energy" type: Real

residual vars
1: fWPH_11_1.nusseltCondenser.p:STATE(1)(min = 611.657 max = 100000000.0 start = fWPH_11_1.flowSTEAM.pstart unit = "Pa" nominal = 1000000.0 stateSelect=StateSelect.prefer )  "Medium pressure" type: Real
2: fWPH_11_1.flowSTEAM.htilde[1]:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.flowSTEAM.hstart[2] unit = "J/kg" nominal = 500000.0 )  "Enthalpy state variables" type: Real [2]
3: fWPH_11_1.flowSTEAM.htilde[2]:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.flowSTEAM.hstart[3] unit = "J/kg" nominal = 500000.0 )  "Enthalpy state variables" type: Real [2]
4: fWPH_11_1.nusseltCondenser.hv:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.nusseltCondenser.hv_start unit = "J/kg" nominal = 500000.0 stateSelect=StateSelect.prefer )  "Specific vapor enthalpy" type: Real
5: $DER.fWPH_11_1.nusseltCondenser.Vv:DUMMY_DER(fixed = false )  "Volume occupied by vapour" type: Real
6: fWPH_11_1.nusseltCondenser.hl:STATE(1)(min = -10000000000.0 max = 10000000000.0 start = fWPH_11_1.nusseltCondenser.hl_start unit = "J/kg" nominal = 500000.0 stateSelect=StateSelect.prefer )  "Specific liquid enthalpy" type: Real

internal equation
1/1 (1): $DER.fWPH_11_1.flowSTEAM.p = der(fWPH_11_1.nusseltCondenser.p)   [dynamic |0|0|0|0|]
2/2 (1): fWPH_11_1.flowSTEAM.dMdt[2] = 4.2 * fWPH_11_1.flowSTEAM.A * (fWPH_11_1.flowSTEAM.drbdh1[2] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.drbdh2[2] * der(fWPH_11_1.flowSTEAM.htilde[2]) + fWPH_11_1.flowSTEAM.drbdp[2] * $DER.fWPH_11_1.flowSTEAM.p)   [dynamic |0|0|0|0|]
3/3 (1): fWPH_11_1.flowSTEAM.dMdt[1] = 4.2 * fWPH_11_1.flowSTEAM.A * (fWPH_11_1.flowSTEAM.drbdh2[1] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.drbdp[1] * $DER.fWPH_11_1.flowSTEAM.p)   [dynamic |0|0|0|0|]
4/4 (1): fWPH_11_1.flowSTEAM.wbar[2] = fWPH_11_1.flowSTEAM.heatTransfer.w[3] + (-0.5) * fWPH_11_1.flowSTEAM.dMdt[2] - fWPH_11_1.flowSTEAM.dMdt[1]   [dynamic |0|0|0|0|]
5/5 (1): fWPH_11_1.flowSTEAM.dMdt[1] + fWPH_11_1.flowSTEAM.dMdt[2] = (fWPH_11_1.STEAMin.m_flow - fWPH_11_1.nusseltCondenser.ws_in) / /*Real*/(fWPH_11_1.flowSTEAM.Nt)   [dynamic |0|0|0|0|]
6/6 (1): fWPH_11_1.flowSTEAM.wbar[1] = fWPH_11_1.flowSTEAM.heatTransfer.w[3] + (-0.5) * fWPH_11_1.flowSTEAM.dMdt[1]   [dynamic |0|0|0|0|]
7/7 (1): $DER.fWPH_11_1.nusseltCondenser.rho_v = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hv, $cse46, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hv))   [unknown |0|0|0|0|]
8/8 (1): $DER.fWPH_11_1.nusseltCondenser.Mv = fWPH_11_1.nusseltCondenser.Vv * $DER.fWPH_11_1.nusseltCondenser.rho_v + $DER.fWPH_11_1.nusseltCondenser.Vv * fWPH_11_1.nusseltCondenser.rho_v   [dynamic |0|0|0|0|]
9/9 (1): $DER.fWPH_11_1.nusseltCondenser.Ev = fWPH_11_1.nusseltCondenser.Mv * (der(fWPH_11_1.nusseltCondenser.hv) + (fWPH_11_1.nusseltCondenser.p * $DER.fWPH_11_1.nusseltCondenser.rho_v - der(fWPH_11_1.nusseltCondenser.p) * fWPH_11_1.nusseltCondenser.rho_v) / fWPH_11_1.nusseltCondenser.rho_v ^ 2.0) + $DER.fWPH_11_1.nusseltCondenser.Mv * (fWPH_11_1.nusseltCondenser.hv - fWPH_11_1.nusseltCondenser.p / fWPH_11_1.nusseltCondenser.rho_v)   [dynamic |0|0|0|0|]
10/10 (1): $DER.fWPH_11_1.nusseltCondenser.rho_l = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hl, $cse47, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hl))   [unknown |0|0|0|0|]
11/11 (1): $DER.fWPH_11_1.nusseltCondenser.Ml = fWPH_11_1.nusseltCondenser.Vl * $DER.fWPH_11_1.nusseltCondenser.rho_l - $DER.fWPH_11_1.nusseltCondenser.Vv * fWPH_11_1.nusseltCondenser.rho_l   [dynamic |0|0|0|0|]
12/12 (1): $DER.fWPH_11_1.nusseltCondenser.El = fWPH_11_1.nusseltCondenser.Ml * (der(fWPH_11_1.nusseltCondenser.hl) + (fWPH_11_1.nusseltCondenser.p * $DER.fWPH_11_1.nusseltCondenser.rho_l - der(fWPH_11_1.nusseltCondenser.p) * fWPH_11_1.nusseltCondenser.rho_l) / fWPH_11_1.nusseltCondenser.rho_l ^ 2.0) + $DER.fWPH_11_1.nusseltCondenser.Ml * (fWPH_11_1.nusseltCondenser.hl - fWPH_11_1.nusseltCondenser.p / fWPH_11_1.nusseltCondenser.rho_l)   [dynamic |0|0|0|0|]

residual equations
1/1 (1): $DER.fWPH_11_1.nusseltCondenser.Ml = fWPH_11_1.nusseltCondenser.ws_in * (1.0 - fWPH_11_1.nusseltCondenser.x_s) + DRAINin.w0 * (1.0 - fWPH_11_1.nusseltCondenser.x_d) + fWPH_11_1.nusseltCondenser.w_cond + (-fWPH_11_1.nusseltCondenser.wd_out) - fWPH_11_1.nusseltCondenser.w_ev   [dynamic |0|0|0|0|]
2/2 (1): $DER.fWPH_11_1.nusseltCondenser.Mv = fWPH_11_1.nusseltCondenser.ws_in * fWPH_11_1.nusseltCondenser.x_s + DRAINin.w0 * fWPH_11_1.nusseltCondenser.x_d + fWPH_11_1.nusseltCondenser.w_ev - fWPH_11_1.nusseltCondenser.w_cond   [dynamic |0|0|0|0|]
3/3 (1): $DER.fWPH_11_1.nusseltCondenser.Ev = fWPH_11_1.nusseltCondenser.ws_in * fWPH_11_1.nusseltCondenser.x_s * $cse53 + DRAINin.w0 * fWPH_11_1.nusseltCondenser.x_d * $cse54 + fWPH_11_1.nusseltCondenser.w_ev * fWPH_11_1.nusseltCondenser.hv_sat + (-fWPH_11_1.nusseltCondenser.Q_cond) - fWPH_11_1.nusseltCondenser.Q_flux   [dynamic |0|0|0|0|]
4/4 (1): 4.2 * fWPH_11_1.flowSTEAM.A * fWPH_11_1.flowSTEAM.rhobar[2] * der(fWPH_11_1.flowSTEAM.htilde[2]) + fWPH_11_1.flowSTEAM.wbar[2] * (fWPH_11_1.flowSTEAM.htilde[2] - fWPH_11_1.flowSTEAM.htilde[1]) + 4.2 * (-fWPH_11_1.flowSTEAM.A) * $DER.fWPH_11_1.flowSTEAM.p = fWPH_11_1.flowSTEAM.Q_single[2]   [dynamic |0|0|0|0|]
5/5 (1): $DER.fWPH_11_1.nusseltCondenser.El = fWPH_11_1.nusseltCondenser.ws_in * (1.0 - fWPH_11_1.nusseltCondenser.x_s) * $cse51 + DRAINin.w0 * (1.0 - fWPH_11_1.nusseltCondenser.x_d) * $cse52 + fWPH_11_1.nusseltCondenser.Q_flux + (-fWPH_11_1.nusseltCondenser.wd_out) * fWPH_11_1.nusseltCondenser.hl - fWPH_11_1.nusseltCondenser.w_ev * fWPH_11_1.nusseltCondenser.hv_sat   [dynamic |0|0|0|0|]
6/6 (1): 4.2 * fWPH_11_1.flowSTEAM.A * fWPH_11_1.flowSTEAM.rho[2] * der(fWPH_11_1.flowSTEAM.htilde[1]) + fWPH_11_1.flowSTEAM.wbar[1] * (fWPH_11_1.flowSTEAM.htilde[1] - STEAMin.h) + 4.2 * (-fWPH_11_1.flowSTEAM.A) * $DER.fWPH_11_1.flowSTEAM.p = fWPH_11_1.flowSTEAM.Q_single[1]   [dynamic |0|0|0|0|]

Unfortunately it is still considered nonlinear and taking approx 10s in total. Looking at the subset of the residual equation related adjacency matrix entries i don't know why it is considered nonlinear:

13:(8,solved) (14,variable(true)) 
14:(11,solved) (14,variable(true)) 
15:(10,solved) (14,variable(true)) 
16:(4,variable(true)) (15,variable(true)) (-4,variable(true)) (-5,variable(true)) (18,param(true)) 
17:(7,solved) (14,variable(true)) (-1,variable(true)) 
18:(5,variable(true)) (13,variable(true)) (-5,variable(true)) (18,param(true))

Maybe there is something basic wrong here and we consider too many systems to be nonlinear, even if they are not.

follow-up: 15 comment:8 by Karim Adbdelhak, 5 years ago

I was able to change the detection of nonlinearity for torn systems such that this is not considered a nonlinear system anymore in PR670. It could still be quite buggy since i was very strict with my algorithms, it quite possibly has to be tweaked to add all exceptions to my rule.

Once again: Unfortunately this did not change much.
BEFORE:

    timeTotal = 9.488673944        

Notification: Torn system details for strict tearing set:
 * Linear torn systems: 14 {(3,88.9%) 4,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 2,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 2
,(1,100.0%) 2,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1}
 * Non-linear torn systems: 7 {1 8,1 8,7 11,1 2,1 1,1 1,1 2}

AFTER:

    timeTotal = 9.25275791
                                              
Notification: Torn system details for strict tearing set:                                                          
 * Linear torn systems: 21 {(3,88.9%) 4,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 2,(1,100.0%) 8,(1,100.0%) 2,(3,88.9%) 4
,(1,100.0%) 2,(1,100.0%) 8,(1,100.0%) 2,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 1,(1,100.0%) 1,(7,46.9%) 11,(1,100.0%)
2,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 2}
 * Non-linear torn systems: 0

Are there supposed to be any nonlinear systems here? If so i might have skipped them. Also this is still with the *wrong* tearing set that celliers algorithm produces. Maybe i can get better results combining it with omcTearing.

follow-up: 10 comment:9 by Francesco Casella, 5 years ago

@Karim, I'll check tomorrow with Dymola, now I can't use it remotely.

I'm a bit confused by what I read in this ticket. Accoring to my understanding, a certain system with certain unknowns x is linear if the Jacobian of the residual functions w.r.t. x does not depend on x. How is it possible that you have different interpretations and different results of such a basic definition?

in reply to: 9 comment:10 by Karim Adbdelhak, 5 years ago

Replying to casella:

@Karim, I'll check tomorrow with Dymola, now I can't use it remotely.

I'm a bit confused by what I read in this ticket. Accoring to my understanding, a certain system with certain unknowns x is linear if the Jacobian of the residual functions w.r.t. x does not depend on x. How is it possible that you have different interpretations and different results of such a basic definition?

The jacobian of the full system is generated to create that solvability information. Since we tear the system we only have to look at a subset of said jacobian, i think we don't do that and look at the full jacobian instead. But i need to do some more analysis, i just created kind of a "hack function" to see if it would change anything. I will try it the correct way tomorrow.

comment:11 by Vitalij Ruge, 5 years ago

Maybe, #3511 is intressting for the issue.

comment:12 by Francesco Casella, 5 years ago

For the record, I understood during today's webmeeting that the "jacobian of the full system" mentioned in comment:10 is computed by some vestigial code in the backend which is now only used for this purpose, but is otherwise no longer used or maintained. @Karim will look into it.

comment:13 by Francesco Casella, 5 years ago

See also #5795, which could be related because it causes strongly connected components much larger than real to be found.

comment:14 by Francesco Casella, 5 years ago

See also #5727 for the original, larger test case.

in reply to: 8 comment:15 by Francesco Casella, 5 years ago

Replying to Karim.Abdelhak:

I was able to change the detection of nonlinearity for torn systems such that this is not considered a nonlinear system anymore in PR670. It could still be quite buggy since i was very strict with my algorithms, it quite possibly has to be tweaked to add all exceptions to my rule.

Once again: Unfortunately this did not change much.
BEFORE:
    timeTotal = 9.488673944        

Notification: Torn system details for strict tearing set:
 * Linear torn systems: 14 {(3,88.9%) 4,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 2,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 2
,(1,100.0%) 2,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1}
 * Non-linear torn systems: 7 {1 8,1 8,7 11,1 2,1 1,1 1,1 2}
AFTER:
    timeTotal = 9.25275791
                                              
Notification: Torn system details for strict tearing set:                                                          
 * Linear torn systems: 21 {(3,88.9%) 4,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 2,(1,100.0%) 8,(1,100.0%) 2,(3,88.9%) 4
,(1,100.0%) 2,(1,100.0%) 8,(1,100.0%) 2,(1,100.0%) 2,(3,88.9%) 4,(1,100.0%) 1,(1,100.0%) 1,(7,46.9%) 11,(1,100.0%)
2,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 2}
 * Non-linear torn systems: 0
Are there supposed to be any nonlinear systems here?

Yes, Dymola finds seven nonlinear systems which are torn to have each a single tearing variable:

fWPH_11_1.flowSTEAM.w
fWPH_11_1.metalTubeFV_cond.ext.T[1]
fWPH_11_1.metalTubeFV_cond.ext.T[2]
fWPH_11_1.flowDRAIN.w
fWPH_11_1.flowFW_cond.infl.m_flow
fWPH_11_1.flowFW_des.infl.m_flow
fWPH_11_1.flowFW_des.wout

If so i might have skipped them.

I guess so. This term shows up in one equation

sqrt(fWPH_11_1.flowFW_des.wout*fWPH_11_1.flowFW_des.wout
      +fWPH_11_1.flowFW_des.wnom*fWPH_11_1.flowFW_des.wnf*(fWPH_11_1.flowFW_des.wnom
      *fWPH_11_1.flowFW_des.wnf)/fWPH_11_1.flowFW_des.Nt/fWPH_11_1.flowFW_des.Nt))

and is clearly nonlinear in fWPH_11_1.flowFW_des.wout

Also this is still with the *wrong* tearing set that celliers algorithm produces. Maybe i can get better results combining it with omcTearing.

No idea.

Last edited 5 years ago by Francesco Casella (previous) (diff)

comment:16 by matthis.thorade@…, 5 years ago

Subscribing

comment:17 by Francesco Casella, 5 years ago

Cc:	matthis.thorade@… added

comment:18 by Francesco Casella, 5 years ago

@Karim reports some progress here.

It involves differentiated functions with dummy inputs that are generating bogus couplings

Version 0, edited 5 years ago by Francesco Casella (next)

in reply to: 5 ; follow-up: 20 comment:19 by Karim Adbdelhak, 5 years ago

It seems to origin from the function Modelica.Media.Water.IF97_Utilities.rho_props_ph (which is used in Modelica.Media.Water.IF97_Utilities.rho_ph) and has a pre defined derivative function Modelica.Media.Water.IF97_Utilities.rho_ph_der.

If we look at the very first dump of comment 5 (between the first main header and the second main header surrounded by ### lines) we can see a loop that has been falsely detected as nonlinear. Linearity analysis is made on the base of the (old) jacobian structure, so i looked into that and realized that the module failed to create a symbolic jacobian (therefore worst case -> nonlinear is assumed and numeric jacobian should be used). I narrowed it down to a differentiation issue that reports that the fifth internal equation equation mentioned in the dump from comment 5

5/5 (1): $DER.fWPH_11_1.nusseltCondenser.rho_l = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hl, $cse47, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hl))   [unknown |0|0|0|0|]

could not be differentiated w.r.t. variable fWPH_11_1.nusseltCondenser.p. The function actually does not use this first input at all, so the derivative should just be zero. Note that this dummy input is necessary to correctly map the inputs from the original function to the derivative function when creating the function call.

Since this should be useful information for many other modules to know which of the inputs are actually used (e.g. adjacency matrix, differentiation, tearing) i created ticket #5874 to collect that information in the frontend and have it right from the get go in the backend.

in reply to: 19 ; follow-up: 21 comment:20 by Francesco Casella, 5 years ago

Replying to Karim.Abdelhak:

I narrowed it down to a differentiation issue that reports that the fifth internal equation equation mentioned in the dump from comment 5
5/5 (1): $DER.fWPH_11_1.nusseltCondenser.rho_l = Modelica.Media.Water.IF97_Utilities.rho_ph_der(fWPH_11_1.nusseltCondenser.p, fWPH_11_1.nusseltCondenser.hl, $cse47, der(fWPH_11_1.nusseltCondenser.p), der(fWPH_11_1.nusseltCondenser.hl))   [unknown |0|0|0|0|]
could not be differentiated w.r.t. variable fWPH_11_1.nusseltCondenser.p. The function actually does not use this first input at all, so the derivative should just be zero.

Wait, are you really sure about that? From what I understand, computing the derivative of the liquid density rho_l does indeed require to know the pressure value fWPH_11_1.nusseltCondenser.p.

From what do you infer that the function is not using the first input at all? As I understand, that shouldn't be the case.

in reply to: 20 comment:21 by Karim Adbdelhak, 5 years ago

Replying to casella:

From what do you infer that the function is not using the first input at all? As I understand, that shouldn't be the case.

This is the function body, taken straight from the MSL:

function rho_ph_der "Derivative function of rho_ph"
  extends Modelica.Icons.Function;
  input SI.Pressure p "Pressure";
  input SI.SpecificEnthalpy h "Specific enthalpy";
  input Common.IF97BaseTwoPhase properties "Auxiliary record";
  input Real p_der "Derivative of pressure";
  input Real h_der "Derivative of specific enthalpy";
  output Real rho_der "Derivative of density";
algorithm
  if (properties.region == 4) then
    rho_der := (properties.rho*(properties.rho*properties.cv/properties.dpT + 1.0)/(properties.dpT*properties.T))*p_der
       + (-properties.rho*properties.rho/(properties.dpT*properties.T))*h_der;
  elseif (properties.region == 3) then
    rho_der := ((properties.rho*(properties.cv*properties.rho + properties.pt))/(properties.rho*properties.rho*properties.pd*
      properties.cv + properties.T*properties.pt*properties.pt))*p_der + (-properties.rho*properties.rho*properties.pt/(properties.rho
      *properties.rho*properties.pd*properties.cv + properties.T*properties.pt*properties.pt))*h_der;
  else
    //regions 1,2,5
    rho_der := (-properties.rho*properties.rho*(properties.vp*properties.cp - properties.vt/properties.rho + properties.T*properties.vt
      *properties.vt)/properties.cp)*p_der + (-properties.rho*properties.rho*properties.vt/(properties.cp))*h_der;
  end if;
end rho_ph_der;

As far as i can see the first input input SI.Pressure p "Pressure"; is not used at all.

Last edited 5 years ago by Karim Adbdelhak (previous) (diff)

follow-up: 23 comment:22 by Francesco Casella, 5 years ago

This is definitely non-trivial. Hubertus designed this very clever mechanism to avoid unecessary repeated calls of the basic heavyweight function waterProps_ph, but it's a nightmare to follow.

This function computes the derivative of the density with respect to pressure and enthalpy, the latter in most cases being strongly related to temperature. So it considers both compressibility and thermal expansion effects.

The compressibility of some media (e.g. gases) obviously depend on the pressure: for the same rate of change of enthalpy or of pressure, a high-pressure gas (which is very dense and thus has large properties.rho) has a much higher density derivative than a low-pressure gas (which has small properties.rho).

Why aren't you taking that indirect dependency into account, when considering the structural dependencies?

in reply to: 22 ; follow-up: 24 comment:23 by Karim Adbdelhak, 5 years ago

Replying to casella:

This is definitely non-trivial. Hubertus designed this very clever mechanism to avoid unecessary repeated calls of the basic heavyweight function waterProps_ph, but it's a nightmare to follow.

This function computes the derivative of the density with respect to pressure and enthalpy, the latter in most cases being strongly related to temperature. So it considers both compressibility and thermal expansion effects.

The compressibility of some media (e.g. gases) obviously depend on the pressure: for the same rate of change of enthalpy or of pressure, a high-pressure gas (which is very dense and thus has large properties.rho) has a much higher density derivative than a low-pressure gas (which has small properties.rho).

I don't think i can follow you here but as far is understand this the function does not directly, but indirectly depend in the pressure.

Why aren't you taking that indirect dependency into account, when considering the structural dependencies?

We do, that is what matching and sorting does. But this function does not contain the variable itself and that has to be considered at many places. E.g. it can never be (implicitely) solved for it even though the inputs depend on it.

But you are right, for the jacobian it also matters.

in reply to: 23 comment:24 by Francesco Casella, 5 years ago

Replying to Karim.Abdelhak:

I don't think i can follow you here but as far is understand this the function does not directly, but indirectly depend in the pressure.

Yes. You have f(x,y) and g(w,z), this is like calling f(x,g(w,z)), which of course depends on x, w and z.

Hence, the total derivative is df/dx*der(x) + df/dy*(dg/dw*der(w) + dg/dz*der(z))

But you are right, for the jacobian it also matters.

That goes without saying.

follow-up: 26 comment:25 by Karim Adbdelhak, 5 years ago

The changes i made produced an error on Modelica.Thermal.FluidHeatFlow.Examples.TwoMass.

I realized that following loop has been considered nonlinear by the old implementation and is now identified to be linear:

Variables:
1: pipe1.flowPort_b.H_flow:VARIABLE(flow=true unit = "W" )  type: Real                                             
2: pipe3.flowPort_a.H_flow:VARIABLE(flow=true unit = "W" )  type: Real                                
3: pipe2.flowPort_b.H_flow:VARIABLE(flow=true unit = "W" )  type: Real                      
4: pipe1.flowPort_b.h:VARIABLE(flow=false unit = "J/kg" )  type: Real                                              

Equations:
1/1 (1): -pipe3.flowPort_a.H_flow = semiLinear(pipe3.flowPort_b.m_flow, pipe3.h, pipe1.flowPort_b.h)   [dynamic |0|
0|0|0|]
2/2 (1): -pipe2.flowPort_b.H_flow = semiLinear(pipe2.flowPort_a.m_flow, pipe2.h, pipe1.flowPort_b.h)   [dynamic |0|
0|0|0|]                       
3/3 (1): pipe1.flowPort_b.H_flow + pipe2.flowPort_b.H_flow + pipe3.flowPort_a.H_flow = 0.0   [dynamic |0|0|0|0|]
4/4 (1): -pipe1.flowPort_b.H_flow = semiLinear(pipe1.flowPort_a.m_flow, pipe1.h, pipe1.flowPort_b.h)   [dynamic |0|
0|0|0|]

I replaced the old jacobian structure with the new one, now it seems to work, i hope it fixes a lot of other errors too.

Just making sure i am not making an error here. This system is linear right? I do not see any reason why it should not be. It is ridiculous to me that it has not been handled as such before.

Before i can push this i need to do a lot of cleaning up work because this should cause a lot of loops to be considered as linear which has to be checked. Also i want to remove the old jacobian structure entirely.

in reply to: 25 comment:26 by Francesco Casella, 5 years ago

Replying to Karim.Abdelhak:

The changes i made produced an error on Modelica.Thermal.FluidHeatFlow.Examples.TwoMass.

I realized that following loop has been considered nonlinear by the old implementation and is now identified to be linear:

If the m_flow variables are known, of course it is :)

I replaced the old jacobian structure with the new one, now it seems to work, i hope it fixes a lot of other errors too.

Wonderful!

Just making sure i am not making an error here. This system is linear right? I do not see any reason why it should not be. It is ridiculous to me that it has not been handled as such before.

I guess semilinear() was identified as nonlinear no matter what. Of course it is only nonlinear w.r.t. its first argument.

Before i can push this i need to do a lot of cleaning up work because this should cause a lot of loops to be considered as linear which has to be checked. Also i want to remove the old jacobian structure entirely.

Sounds good. I hope we'll have a significant improvement in robustness and speed once this is done.

Thanks for the hard work!

comment:27 by Francesco Casella, 4 years ago

Milestone:	1.16.0 → 1.17.0

Retargeted to 1.17.0 after 1.16.0 release

comment:28 by Francesco Casella, 4 years ago

Milestone:	1.17.0 → 1.18.0

Retargeted to 1.18.0 because of 1.17.0 timed release.

comment:29 by Francesco Casella, 3 years ago

Milestone:	1.18.0

Ticket retargeted after milestone closed

Note: See TracTickets for help on using tickets.

Download in other formats: