Model lag for inductor capacitor circuit

In https://github.com/phetsims/circuit-construction-kit-common/issues/772#issuecomment-965690073 @stemilymill said:

For https://github.com/phetsims/qa/issues/736

I'm not sure if this is helpful but it seems related:

Several times I was able to freeze the sim in the manner shown in this video. You can still see the current still flowing here but the sim is not responsive to user input and eventually completely crashes. This was taken with a very low frame rate, and I sped up the beginning (building the circuit) so the choppiness is not actual performance. I can work on getting a higher framerate recording.

https://user-images.githubusercontent.com/85511101/141181334-c6f1c341-487e-4651-b6fb-3dc4f92b4bae.mp4

I found that increasing ?minDT from the default value of ?minDT=1E-3 to ?minDT=1E-2 corrected this problem. So it does seem like another occurrence of poor matrix conditioning in an LC circuit due to low dt values. However, we tried to turn the minDT as low as we could in order to get accurate physics. We recently turned up minDT because when it was too low we suffered from performance problems. I'm at a good point to check in with @arouinfar and @ariel-phet about what's possible.

@arouinfar and I discussed this today, we were able to reproduce the problem and see it disappear when dt was increased to 1E-2. However, changing the dt to 1E-2 causes some values to go past our previously determined error threshold of 0.05 in the unit tests. We discussed combining adjacent series inductors or capacitors into a single composite for purposes of the circuit solution--we were hopeful this may work because in the circuit above, we did not see a problem when we removed one inductor from the loop. However, this is a complex strategy that could take a week or so and would introduce questions like how to backpropagate solutions from the composite model back to individual circuit elements.

We also would like to see how the unit tests are failing when dt goes to 1E-2. If it is failing from a trailing or leading temporal problem, maybe it will be ok, since the user cannot start a clock in CCK within a fine resolution anyways. So I'll test that next and @arouinfar and I will discuss again shortly.

According to http://dev.hypertriton.com/edacious/trunk/doc/lec.pdf we can follow the Norton model for the inductor companion model to have all of the dt in the numerator. But we will have to add current sources. It's a current source in parallel with a resistor, with:

Req = dt/2/L
Ieq = i0 + Req * v0

This seems like the right thing to investigate next, since we would like to be able to turn down the dt without creating a poorly conditioned problem.

I implemented the norton inductor trapezoidal rule. Not sure about all the signs, but unit tests are failing.

```diff Index: js/model/analysis/LTACircuit.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/LTACircuit.ts b/js/model/analysis/LTACircuit.ts --- a/js/model/analysis/LTACircuit.ts (revision 4dd54012559e66542092a7ac8b41239c5428be08) +++ b/js/model/analysis/LTACircuit.ts (date 1637039447790) @@ -55,6 +55,7 @@ const companionBatteries: MNABattery[] = []; const companionResistors: MNAResistor[] = []; + const companionCurrents: MNACurrent[] = []; const currentCompanions: { element: CoreModel, getValueForSolution: ( solution: MNASolution ) => number }[] = []; // Node indices that have been used @@ -124,35 +125,49 @@ // See najm page 279 and Pillage page 86 this.ltaInductors.forEach( ltaInductor => { - // In series - const newNode = 'syntheticNode' + ( syntheticNodeIndex++ ); - const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); + // Req = dt/2/L + // Ieq = i0 + Req * v0 + const companionResistance = dt / 2 / ltaInductor.inductance; + const companionCurrent = ltaInductor.current + companionResistance * ltaInductor.voltage; - const companionResistance = 2 * ltaInductor.inductance / dt; - const companionVoltage = -ltaInductor.voltage - companionResistance * ltaInductor.current; - - const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); - const resistor = new MNAResistor( newNode, newNode2, companionResistance ); - const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); - companionBatteries.push( battery ); + const resistor = new MNAResistor( ltaInductor.node0, ltaInductor.node1, companionResistance ); companionResistors.push( resistor ); - companionResistors.push( addedResistor ); + companionCurrents.push( new MNACurrent( ltaInductor.node0, ltaInductor.node1, -companionCurrent ) ); + ltaInductor.inductorVoltageNode1 = ltaInductor.node1; - ltaInductor.inductorVoltageNode1 = newNode2; - - // we need to be able to get the current for this component - // in series, so current is same through both companion components currentCompanions.push( { element: ltaInductor, - // TODO: This sign looks very wrong https://github.com/phetsims/circuit-construction-kit-common/issues/758 - getValueForSolution: ( solution: MNASolution ) => -solution.getCurrentForResistor( resistor ) + getValueForSolution: ( solution: MNASolution ) => solution.getCurrentForResistor( resistor ) + companionCurrent } ); + + // In series + // const newNode = 'syntheticNode' + ( syntheticNodeIndex++ ); + // const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); + // + // const companionResistance = 2 * ltaInductor.inductance / dt; + // const companionVoltage = -ltaInductor.voltage - companionResistance * ltaInductor.current; + // + // const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); + // const resistor = new MNAResistor( newNode, newNode2, companionResistance ); + // const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); + // companionBatteries.push( battery ); + // companionResistors.push( resistor ); + // companionResistors.push( addedResistor ); + // + // ltaInductor.inductorVoltageNode1 = newNode2; + // + // // we need to be able to get the current for this component + // // in series, so current is same through both companion components + // currentCompanions.push( { + // element: ltaInductor, + // // TODO: This sign looks very wrong https://github.com/phetsims/circuit-construction-kit-common/issues/758 + // getValueForSolution: ( solution: MNASolution ) => -solution.getCurrentForResistor( resistor ) + // } ); } ); const newBatteryList = companionBatteries; const newResistorList = [ ...this.ltaResistors, ...companionResistors ]; - const newCurrentList: MNACurrent[] = []; // Placeholder for if we add other circuit elements in the future - + const newCurrentList = companionCurrents; const mnaCircuit = new MNACircuit( newBatteryList, newResistorList, newCurrentList ); const mnaSolution = mnaCircuit.solve(); @@ -235,7 +250,7 @@ ltaInductor.node0, ltaInductor.node1, solution.getVoltage( ltaInductor.node0, ltaInductor.inductorVoltageNode1! ), - solution.getCurrentForCompanion( ltaInductor ), + solution.getCurrentForCompanion( ltaInductor ), // TODO: Should this be the whole current or just the part through the ltaInductor.inductance ); } ); ```

To test:

Did the unit tests have the wrong signs?
I added a test from khan academy with a current source, and have some sign questions there as well.

I investigated all combinations of signs in this patch, none gave expected behavior (even ignoring unit tests signs). Next we should go back to the standard MNA and check all the signs in case that has an impact here.

```diff Index: js/model/analysis/LTACircuitTests.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/LTACircuitTests.ts b/js/model/analysis/LTACircuitTests.ts --- a/js/model/analysis/LTACircuitTests.ts (revision 4dd54012559e66542092a7ac8b41239c5428be08) +++ b/js/model/analysis/LTACircuitTests.ts (date 1637040931058) @@ -24,64 +24,67 @@ let worstError = 0; // @ts-ignore -// window.string = `{v} {r} {c} {t} {desiredVoltageAtTPlusDT} {voltage} {error} -// `; +window.capacitorString = `{v} {r} {c} {t} {desiredVoltageAtTPlusDT} {voltage} {error} +`; +// @ts-ignore +window.inductorString = `{v} {r} {L} {t} {expectedCurrent} {actualCurrent} {error} +`; // Check the values coming from an RC circuit (may be an equivalent capacitance) -const iterateCapacitor = ( circuit: LTACircuit, resistor: MNAResistor, v: number, r: number, c: number, assert: Assert ) => { - for ( let i = 0; i < ITERATIONS; i++ ) { - const t = i * dt; - - const companionSolution = circuit.solveItWithSubdivisions( dt ); - const voltage = companionSolution!.getVoltage( resistor.nodeId0, resistor.nodeId1 ); - const desiredVoltageAtTPlusDT = -v * Math.exp( -( t + dt ) / r / c ); - const error = Math.abs( voltage - desiredVoltageAtTPlusDT ); - - // console.log( error ); - if ( error > worstError ) { - worstError = error; - // console.log( 'new worst error: ' + worstError ); - } - - // @ts-ignore - // window.string = window.string + `${v} ${r} ${c} ${t} ${desiredVoltageAtTPlusDT} ${voltage} ${error} -// `; - assert.ok( error < errorThreshold ); // sample run indicates largest error is 1.5328E-7 - circuit = circuit.updateWithSubdivisions( dt ); - } -}; - -const testVRCCircuit = ( v: number, r: number, c: number, assert: Assert ) => { - const resistor = new MNAResistor( '1', '2', r ); - const battery = new LTAResistiveBattery( id++, '0', '1', v, 0 ); - const capacitor = new LTACapacitor( id++, '2', '0', 0.0, v / r, c ); - const circuit = new LTACircuit( [ resistor ], [ battery ], [ capacitor ], [] ); - iterateCapacitor( circuit, resistor, v, r, c, assert ); -}; - -const testVRCCircuitSeriesCapacitors = ( v: number, r: number, c1: number, c2: number, assert: Assert ) => { - - const ceq = 1 / ( 1 / c1 + 1 / c2 ); - const resistor = new MNAResistor( '1', '2', r ); - const battery = new LTAResistiveBattery( id++, '0', '1', v, 0 ); - const capacitor1 = new LTACapacitor( id++, '2', '3', 0.0, v / r, c1 ); - const capacitor2 = new LTACapacitor( id++, '3', '0', 0.0, v / r, c2 ); - - const circuit = new LTACircuit( [ resistor ], [ battery ], [ capacitor1, capacitor2 ], [] ); - - iterateCapacitor( circuit, resistor, v, r, ceq, assert ); -}; - -const testVRCCircuitParallelCapacitors = ( v: number, r: number, c1: number, c2: number, assert: Assert ) => { - - const ceq = c1 + c2; - const resistor = new MNAResistor( '1', '2', r ); - const battery = new LTAResistiveBattery( id++, '0', '1', v, 0 ); - const capacitor1 = new LTACapacitor( id++, '2', '0', 0.0, v / r, c1 ); - const capacitor2 = new LTACapacitor( id++, '2', '0', 0.0, v / r, c2 ); - - const circuit = new LTACircuit( [ resistor ], [ battery ], [ capacitor1, capacitor2 ], [] ); - iterateCapacitor( circuit, resistor, v, r, ceq, assert ); -}; +// const iterateCapacitor = ( circuit: LTACircuit, resistor: MNAResistor, v: number, r: number, c: number, assert: Assert ) => { +// for ( let i = 0; i < ITERATIONS; i++ ) { +// const t = i * dt; +// +// const companionSolution = circuit.solveItWithSubdivisions( dt ); +// const voltage = companionSolution!.getVoltage( resistor.nodeId0, resistor.nodeId1 ); +// const desiredVoltageAtTPlusDT = -v * Math.exp( -( t + dt ) / r / c ); +// const error = Math.abs( voltage - desiredVoltageAtTPlusDT ); +// +// // console.log( error ); +// if ( error > worstError ) { +// worstError = error; +// // console.log( 'new worst error: ' + worstError ); +// } +// +// // @ts-ignore +// // window.capacitorString = window.capacitorString + `${v} ${r} ${c} ${t} ${desiredVoltageAtTPlusDT} ${voltage} ${error} +// // `; +// assert.ok( error < errorThreshold ); // sample run indicates largest error is 1.5328E-7 +// circuit = circuit.updateWithSubdivisions( dt ); +// } +// }; +// +// const testVRCCircuit = ( v: number, r: number, c: number, assert: Assert ) => { +// const resistor = new MNAResistor( '1', '2', r ); +// const battery = new LTAResistiveBattery( id++, '0', '1', v, 0 ); +// const capacitor = new LTACapacitor( id++, '2', '0', 0.0, v / r, c ); +// const circuit = new LTACircuit( [ resistor ], [ battery ], [ capacitor ], [] ); +// iterateCapacitor( circuit, resistor, v, r, c, assert ); +// }; +// +// const testVRCCircuitSeriesCapacitors = ( v: number, r: number, c1: number, c2: number, assert: Assert ) => { +// +// const ceq = 1 / ( 1 / c1 + 1 / c2 ); +// const resistor = new MNAResistor( '1', '2', r ); +// const battery = new LTAResistiveBattery( id++, '0', '1', v, 0 ); +// const capacitor1 = new LTACapacitor( id++, '2', '3', 0.0, v / r, c1 ); +// const capacitor2 = new LTACapacitor( id++, '3', '0', 0.0, v / r, c2 ); +// +// const circuit = new LTACircuit( [ resistor ], [ battery ], [ capacitor1, capacitor2 ], [] ); +// +// iterateCapacitor( circuit, resistor, v, r, ceq, assert ); +// }; +// +// const testVRCCircuitParallelCapacitors = ( v: number, r: number, c1: number, c2: number, assert: Assert ) => { +// +// const ceq = c1 + c2; +// const resistor = new MNAResistor( '1', '2', r ); +// const battery = new LTAResistiveBattery( id++, '0', '1', v, 0 ); +// const capacitor1 = new LTACapacitor( id++, '2', '0', 0.0, v / r, c1 ); +// const capacitor2 = new LTACapacitor( id++, '2', '0', 0.0, v / r, c2 ); +// +// const circuit = new LTACircuit( [ resistor ], [ battery ], [ capacitor1, capacitor2 ], [] ); +// iterateCapacitor( circuit, resistor, v, r, ceq, assert ); +// }; // This is for comparison with TestTheveninCapacitorRC // QUnit.test( 'testVRC991Eminus2', assert => { @@ -93,35 +96,35 @@ // } // } ); -QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 5 r 10 c 1E minus2', assert => { - testVRCCircuit( 5.0, 10.0, 1.0E-2, assert ); -} ); - -QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 10 r 10 c 1E minus2', assert => { - testVRCCircuit( 10.0, 10.0, 1.0E-2, assert ); -} ); - -QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 3 r 7 c 1Eminus1', assert => { - testVRCCircuit( 3, 7, 1E-1, assert ); -} ); - -QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 3 r 7 c 100', assert => { - testVRCCircuit( 3, 7, 100, assert ); -} ); - -QUnit.test( 'test RC Circuit with series capacitors', assert => { - testVRCCircuitSeriesCapacitors( 3, 7, 10, 10, assert ); - for ( let i = 0; i < 10; i++ ) { - testVRCCircuitSeriesCapacitors( 3, 7, Math.random() * 10, Math.random() * 10, assert ); // eslint-disable-line - } -} ); - -QUnit.test( 'test RC Circuit with parallel capacitors', assert => { - testVRCCircuitParallelCapacitors( 3, 7, 10, 10, assert ); - for ( let i = 0; i < 10; i++ ) { - testVRCCircuitParallelCapacitors( 3, 7, Math.random() * 10, Math.random() * 10, assert ); // eslint-disable-line - } -} ); +// QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 5 r 10 c 1E minus2', assert => { +// testVRCCircuit( 5.0, 10.0, 1.0E-2, assert ); +// } ); +// +// QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 10 r 10 c 1E minus2', assert => { +// testVRCCircuit( 10.0, 10.0, 1.0E-2, assert ); +// } ); +// +// QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 3 r 7 c 1Eminus1', assert => { +// testVRCCircuit( 3, 7, 1E-1, assert ); +// } ); +// +// QUnit.test( 'test RC Circuit should have voltage exponentially decay with T RC for v 3 r 7 c 100', assert => { +// testVRCCircuit( 3, 7, 100, assert ); +// } ); +// +// QUnit.test( 'test RC Circuit with series capacitors', assert => { +// testVRCCircuitSeriesCapacitors( 3, 7, 10, 10, assert ); +// for ( let i = 0; i < 10; i++ ) { +// testVRCCircuitSeriesCapacitors( 3, 7, Math.random() * 10, Math.random() * 10, assert ); // eslint-disable-line +// } +// } ); +// +// QUnit.test( 'test RC Circuit with parallel capacitors', assert => { +// testVRCCircuitParallelCapacitors( 3, 7, 10, 10, assert ); +// for ( let i = 0; i < 10; i++ ) { +// testVRCCircuitParallelCapacitors( 3, 7, Math.random() * 10, Math.random() * 10, assert ); // eslint-disable-line +// } +// } ); const iterateInductor = ( circuit: LTACircuit, resistor: MNAResistor, V: number, R: number, L: number, assert: Assert ) => { for ( let i = 0; i < ITERATIONS; i++ ) { @@ -130,7 +133,13 @@ const current = solution!.getCurrent( resistor )!; const expectedCurrent = V / R * ( 1 - Math.exp( -( t + dt ) * R / L ) );//positive, by definition of MNA.Battery const error = Math.abs( current - expectedCurrent ); - assert.ok( error < errorThreshold ); + + // @ts-ignore + window.inductorString = window.inductorString + `${V} ${R} ${L} ${t} ${expectedCurrent} ${current} ${error} +`; + console.log( dt, expectedCurrent, current ); + // assert.ok( error < errorThreshold ); + assert.ok( true ); circuit = circuit.updateWithSubdivisions( dt ); } }; @@ -142,42 +151,44 @@ const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor ] ); iterateInductor( circuit, resistor, V, R, L, assert ); }; -const testVRLCircuitSeries = ( V: number, R: number, L1: number, L2: number, assert: Assert ) => { - const Leq = L1 + L2; - const resistor = new MNAResistor( '1', '2', R ); - const battery = new LTAResistiveBattery( id++, '0', '1', V, 0 ); - const inductor1 = new LTAInductor( id++, '2', '3', V, 0.0, L1 ); - const inductor2 = new LTAInductor( id++, '3', '0', V, 0.0, L2 ); - const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor1, inductor2 ] ); - iterateInductor( circuit, resistor, V, R, Leq, assert ); -}; -const testVRLCircuitParallel = ( V: number, R: number, L1: number, L2: number, assert: Assert ) => { - const Leq = 1 / ( 1 / L1 + 1 / L2 ); - const resistor = new MNAResistor( '1', '2', R ); - const battery = new LTAResistiveBattery( id++, '0', '1', V, 0 ); - const inductor1 = new LTAInductor( id++, '2', '0', V, 0.0, L1 ); - const inductor2 = new LTAInductor( id++, '2', '0', V, 0.0, L2 ); - const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor1, inductor2 ] ); - iterateInductor( circuit, resistor, V, R, Leq, assert ); -}; - +// const testVRLCircuitSeries = ( V: number, R: number, L1: number, L2: number, assert: Assert ) => { +// const Leq = L1 + L2; +// const resistor = new MNAResistor( '1', '2', R ); +// const battery = new LTAResistiveBattery( id++, '0', '1', V, 0 ); +// const inductor1 = new LTAInductor( id++, '2', '3', V, 0.0, L1 ); +// const inductor2 = new LTAInductor( id++, '3', '0', V, 0.0, L2 ); +// const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor1, inductor2 ] ); +// iterateInductor( circuit, resistor, V, R, Leq, assert ); +// }; +// const testVRLCircuitParallel = ( V: number, R: number, L1: number, L2: number, assert: Assert ) => { +// const Leq = 1 / ( 1 / L1 + 1 / L2 ); +// const resistor = new MNAResistor( '1', '2', R ); +// const battery = new LTAResistiveBattery( id++, '0', '1', V, 0 ); +// const inductor1 = new LTAInductor( id++, '2', '0', V, 0.0, L1 ); +// const inductor2 = new LTAInductor( id++, '2', '0', V, 0.0, L2 ); +// const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor1, inductor2 ] ); +// iterateInductor( circuit, resistor, V, R, Leq, assert ); +// }; +// QUnit.test( 'test_RL_Circuit_should_have_correct_behavior', assert => { testVRLCircuit( 5, 10, 1, assert ); - testVRLCircuit( 3, 11, 2.5, assert ); - testVRLCircuit( 7, 13, 1E4, assert ); - testVRLCircuit( 7, 13, 1E-1, assert ); -} ); - -QUnit.test( 'Series inductors', assert => { - testVRLCircuitSeries( 5, 10, 1, 1, assert ); - for ( let i = 0; i < 10; i++ ) { - testVRLCircuitSeries( 10, 10, 5 * Math.random() + 0.1, 5 * Math.random() + 0.1, assert )// eslint-disable-line - } + // testVRLCircuit( 3, 11, 2.5, assert ); + // testVRLCircuit( 7, 13, 1E4, assert ); + // testVRLCircuit( 7, 13, 1E-1, assert ); } ); -QUnit.test( 'Parallel inductors', assert => { - testVRLCircuitParallel( 5, 10, 1, 1, assert ); - for ( let i = 0; i < 10; i++ ) { - testVRLCircuitParallel( 10, 10, 5 * Math.random() + 1, 5 * Math.random() + 1, assert )// eslint-disable-line - } -} ); \ No newline at end of file +// console.log( window.inductorString ); +// +// QUnit.test( 'Series inductors', assert => { +// testVRLCircuitSeries( 5, 10, 1, 1, assert ); +// for ( let i = 0; i < 10; i++ ) { +// testVRLCircuitSeries( 10, 10, 5 * Math.random() + 0.1, 5 * Math.random() + 0.1, assert )// eslint-disable-line +// } +// } ); +// +// QUnit.test( 'Parallel inductors', assert => { +// testVRLCircuitParallel( 5, 10, 1, 1, assert ); +// for ( let i = 0; i < 10; i++ ) { +// testVRLCircuitParallel( 10, 10, 5 * Math.random() + 1, 5 * Math.random() + 1, assert )// eslint-disable-line +// } +// } ); \ No newline at end of file Index: js/model/analysis/LTACircuit.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/LTACircuit.ts b/js/model/analysis/LTACircuit.ts --- a/js/model/analysis/LTACircuit.ts (revision 4dd54012559e66542092a7ac8b41239c5428be08) +++ b/js/model/analysis/LTACircuit.ts (date 1637041727404) @@ -55,6 +55,7 @@ const companionBatteries: MNABattery[] = []; const companionResistors: MNAResistor[] = []; + const companionCurrents: MNACurrent[] = []; const currentCompanions: { element: CoreModel, getValueForSolution: ( solution: MNASolution ) => number }[] = []; // Node indices that have been used @@ -124,35 +125,55 @@ // See najm page 279 and Pillage page 86 this.ltaInductors.forEach( ltaInductor => { - // In series - const newNode = 'syntheticNode' + ( syntheticNodeIndex++ ); - const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); + // Norton Trapezoidal Inductor, see http://dev.hypertriton.com/edacious/trunk/doc/lec.pdf + // Req = dt/2/L + // Ieq = i0 + Req * v0 - const companionResistance = 2 * ltaInductor.inductance / dt; - const companionVoltage = -ltaInductor.voltage - companionResistance * ltaInductor.current; + const A = 1; + const B = 1; + const C = 1; + const D = 1; + const companionResistance = dt / 2 / ltaInductor.inductance; + const companionCurrent = A * ltaInductor.current + B * companionResistance * ltaInductor.voltage; - const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); - const resistor = new MNAResistor( newNode, newNode2, companionResistance ); - const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); - companionBatteries.push( battery ); + const resistor = new MNAResistor( ltaInductor.node0, ltaInductor.node1, companionResistance ); companionResistors.push( resistor ); - companionResistors.push( addedResistor ); + companionCurrents.push( new MNACurrent( ltaInductor.node0, ltaInductor.node1, companionCurrent ) ); + ltaInductor.inductorVoltageNode1 = ltaInductor.node1; - ltaInductor.inductorVoltageNode1 = newNode2; - - // we need to be able to get the current for this component - // in series, so current is same through both companion components currentCompanions.push( { element: ltaInductor, - // TODO: This sign looks very wrong https://github.com/phetsims/circuit-construction-kit-common/issues/758 - getValueForSolution: ( solution: MNASolution ) => -solution.getCurrentForResistor( resistor ) + getValueForSolution: ( solution: MNASolution ) => C * solution.getCurrentForResistor( resistor ) + D * companionCurrent } ); + + // In series + // const newNode = 'syntheticNode' + ( syntheticNodeIndex++ ); + // const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); + // + // const companionResistance = 2 * ltaInductor.inductance / dt; + // const companionVoltage = -ltaInductor.voltage - companionResistance * ltaInductor.current; + // + // const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); + // const resistor = new MNAResistor( newNode, newNode2, companionResistance ); + // const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); + // companionBatteries.push( battery ); + // companionResistors.push( resistor ); + // companionResistors.push( addedResistor ); + // + // ltaInductor.inductorVoltageNode1 = newNode2; + // + // // we need to be able to get the current for this component + // // in series, so current is same through both companion components + // currentCompanions.push( { + // element: ltaInductor, + // // TODO: This sign looks very wrong https://github.com/phetsims/circuit-construction-kit-common/issues/758 + // getValueForSolution: ( solution: MNASolution ) => -solution.getCurrentForResistor( resistor ) + // } ); } ); const newBatteryList = companionBatteries; const newResistorList = [ ...this.ltaResistors, ...companionResistors ]; - const newCurrentList: MNACurrent[] = []; // Placeholder for if we add other circuit elements in the future - + const newCurrentList = companionCurrents; const mnaCircuit = new MNACircuit( newBatteryList, newResistorList, newCurrentList ); const mnaSolution = mnaCircuit.solve(); @@ -235,7 +256,7 @@ ltaInductor.node0, ltaInductor.node1, solution.getVoltage( ltaInductor.node0, ltaInductor.inductorVoltageNode1! ), - solution.getCurrentForCompanion( ltaInductor ), + solution.getCurrentForCompanion( ltaInductor ), // TODO: Should this be the whole current or just the part through the ltaInductor.inductance ); } ); ```

@arouinfar and @ariel-phet and I worked on this today. We tried out the norton equivalent trapezoidal inductor and could not get it to work. We tried all combinations of signs and saw 3 behaviors:

steady -0.5A
oscillating exponential growth
non-oscillating exponential growth

We also saw that we need to seed the inductors with the DC operating values.

We also tried the forward euler, but could not look up the current through the resistor. We tried backward euler and had the same behavior as trapezoidal.

We discussed testing some of the buggy circumstances in the spice simulation. We tried running an RL circuit, but got the error Error: no data saved for Transient analysis; analysis not run. We also discussed how we might wire up the spice simulator into the sim. This is discussed more in https://github.com/phetsims/circuit-construction-kit-common/issues/228

Page 25 of https://docs.lib.purdue.edu/cgi/viewcontent.cgi?article=1301&context=ecetr also has a treatment of inductor companion models. That page shows that for forward euler, there is no parallel resistor. I tried this:

```ts // Copyright 2019-2021, University of Colorado Boulder /** * There are two parts to solving a dynamic circuit: * 1. Splitting up dynamic components such as capacitors and inductors into their respective linear companion models. * 2. Adjusting the dt so that integration steps are accurate. This is done through TimestepSubdivisions. * * @author Sam Reid (PhET Interactive Simulations) */ import CCKCQueryParameters from '../../CCKCQueryParameters.js'; import LTAStateSet from './LTAStateSet.js'; import circuitConstructionKitCommon from '../../circuitConstructionKitCommon.js'; import MNACircuit from './mna/MNACircuit.js'; import TimestepSubdivisions from './TimestepSubdivisions.js'; import LTASolution from './LTASolution.js'; import LTAState from './LTAState.js'; import LTAInductor from './LTAInductor.js'; import LTACapacitor from './LTACapacitor.js'; import MNASolution from './mna/MNASolution.js'; import LTAResistiveBattery from './LTAResistiveBattery.js'; import MNABattery from './mna/MNABattery.js'; import MNAResistor from './mna/MNAResistor.js'; import MNACurrent from './mna/MNACurrent.js'; import CoreModel from './CoreModel.js'; type DistanceParams = { getCharacteristicArray: () => number[]; }; class LTACircuit { private readonly ltaResistors: MNAResistor[]; private readonly ltaBatteries: LTAResistiveBattery[]; readonly ltaCapacitors: LTACapacitor[]; readonly ltaInductors: LTAInductor[]; constructor( ltaResistors: MNAResistor[], ltaBatteries: LTAResistiveBattery[], ltaCapacitors: LTACapacitor[], ltaInductors: LTAInductor[] ) { // @private this.ltaResistors = ltaResistors; this.ltaBatteries = ltaBatteries; this.ltaCapacitors = ltaCapacitors; this.ltaInductors = ltaInductors; } /** * Solving the companion model is the same as propagating forward in time by dt. * * @param {number} dt * @returns {LTASolution} * @public */ solvePropagate( dt: number ) { const companionBatteries: MNABattery[] = []; const companionResistors: MNAResistor[] = []; const companionCurrents: MNACurrent[] = []; const currentCompanions: { element: CoreModel, getValueForSolution: ( solution: MNASolution ) => number }[] = []; // Node indices that have been used let syntheticNodeIndex = 0; // Each resistive battery is a resistor in series with a battery this.ltaBatteries.forEach( resistiveBatteryAdapter => { const newNode = 'syntheticNode' + ( syntheticNodeIndex++ ); const idealBattery = new MNABattery( resistiveBatteryAdapter.node0, newNode, resistiveBatteryAdapter.voltage ); // final LinearCircuitSolver.Battery const idealResistor = new MNAResistor( newNode, resistiveBatteryAdapter.node1, resistiveBatteryAdapter.resistance ); // LinearCircuitSolver.Resistor companionBatteries.push( idealBattery ); companionResistors.push( idealResistor ); // We need to be able to get the current for this component currentCompanions.push( { element: resistiveBatteryAdapter, getValueForSolution: ( solution: MNASolution ) => { return solution.getSolvedCurrent( idealBattery ); } } ); } ); // Add companion models for capacitor // TRAPEZOIDAL: battery and resistor in series. // We use trapezoidal rather than backward Euler because we do not model current sources and it seems to work well. // See http://circsimproj.blogspot.com/2009/07/companion-models.html // Veq = V + dt*I/2/C; // Req = dt/2/C this.ltaCapacitors.forEach( ltaCapacitor => { assert && assert( dt >= 0, 'dt should be non-negative' ); const newNode1 = 'syntheticNode' + ( syntheticNodeIndex++ ); const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); const companionResistance = dt / 2.0 / ltaCapacitor.capacitance; const resistanceTerm = CCKCQueryParameters.capacitorResistance; // The capacitor is modeled as a battery in series with a resistor. Hence the voltage drop across the capacitor // is equal to the voltage drop across the battery plus the voltage drop across the resistor. // V = Vbattery + Vresistor. We need to solve for the voltage across the battery to use it in the companion // model, so we have Vbattery = V-Vresistor. The magnitude of the voltage drop across the resistor is given by // |V|=|IReq| and sign is unchanged since the conventional current flows from high to low voltage. const companionVoltage = companionResistance * ltaCapacitor.current - ltaCapacitor.voltage; const battery = new MNABattery( ltaCapacitor.node0, newNode1, companionVoltage ); const resistor = new MNAResistor( newNode1, newNode2, companionResistance ); const resistor2 = new MNAResistor( newNode2, ltaCapacitor.node1, resistanceTerm ); companionBatteries.push( battery ); companionResistors.push( resistor ); companionResistors.push( resistor2 ); ltaCapacitor.capacitorVoltageNode1 = newNode2; // We need to be able to get the current for this component. In series, so the current is the same through both. currentCompanions.push( { element: ltaCapacitor, getValueForSolution: ( solution: MNASolution ) => solution.getCurrentForResistor( resistor ) } ); } ); const signs = [ 1, 1, 1, 1 ]; // const signs = [ -1, +1, +1, +1 ]; // const signs = [ +1, -1, +1, +1 ]; // const signs = [ +1, 1, -1, +1 ]; // const signs = [ +1, 1, 1, -1 ]; // const signs = [ -1, -1, 1, 1 ]; // const signs = [ -1, 1, -1, 1 ]; // const signs = [ -1, 1, 1, -1 ]; // const signs = [ 1, -1, -1, 1 ]; // const signs = [ 1, -1, 1, -1 ]; // const signs = [ 1, 1, -1, -1 ]; // const signs = [ 1, -1, -1, -1 ]; // const signs = [ -1, 1, -1, -1 ]; // const signs = [ -1, -1, 1, -1 ]; // const signs = [ -1, -1, -1, 1 ]; // const signs = [ -1, -1, -1, -1 ]; // const A = +1; // const B = +1; // const C = +1; // const D = +1; // See also http://circsimproj.blogspot.com/2009/07/companion-models.html, which reports: // Req = 2L/dt // Veq = -2Li/dt-v // See najm page 279 and Pillage page 86 this.ltaInductors.forEach( ltaInductor => { // debugger; // const newNode1 = 'syntheticNode' + ( syntheticNodeIndex++ ); // const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); // Req = dt/2/L // Ieq = i0 + Req * v0 // const companionResistance = dt / ltaInductor.inductance; const companionCurrent = -dt * ltaInductor.voltage / ltaInductor.inductance + ltaInductor.current; // signs[ 0 ] * ltaInductor.current + signs[ 1 ] * companionResistance * ltaInductor.voltage; // console.log( companionCurrent ); // const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); // const resistor = new MNAResistor( ltaInductor.node0, ltaInductor.node1, companionResistance ); companionCurrents.push( new MNACurrent( ltaInductor.node0, ltaInductor.node1, companionCurrent ) ); // const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); // companionBatteries.push( battery ); // companionResistors.push( resistor ); ltaInductor.inductorVoltageNode1 = ltaInductor.node1; currentCompanions.push( { element: ltaInductor, getValueForSolution: ( solution: MNASolution ) => { return companionCurrent; // const a = signs[ 2 ] * solution.getCurrentForResistor( resistor ); // return a + signs[ 3 ] * companionCurrent; } } ); } ); const newBatteryList = companionBatteries; const newResistorList = [ ...this.ltaResistors, ...companionResistors ]; const newCurrentList = companionCurrents; const mnaCircuit = new MNACircuit( newBatteryList, newResistorList, newCurrentList ); const mnaSolution = mnaCircuit.solve(); return new LTASolution( this, mnaSolution, currentCompanions ); } /** * @param {TimestepSubdivisions} timestepSubdivisions * @param {number} dt * @returns {LTAStateSet} * @public */ solveWithSubdivisions( timestepSubdivisions: TimestepSubdivisions, dt: number ) { const steppable = { update: ( a: { update: ( dt: number ) => LTAState; }, dt: number ) => a.update( dt ), distance: ( a: DistanceParams, b: DistanceParams ) => euclideanDistance( a.getCharacteristicArray(), b.getCharacteristicArray() ) }; // Turning the error threshold too low here can fail the inductor tests in MNATestCase const x = timestepSubdivisions.stepInTimeWithHistory( new LTAState( this, null ), steppable, dt ); return new LTAStateSet( x ); } /** * @param {number} dt * @returns {LTAStateSet} * @private */ solveWithSubdivisions2( dt: number ) { return this.solveWithSubdivisions( new TimestepSubdivisions(), dt ); } /** * @param {number} dt * @returns {LTACircuit} * @private */ updateWithSubdivisions( dt: number ) { return this.solveWithSubdivisions2( dt ).getFinalState().ltaCircuit; } /** * @param {number} dt * @returns {LTASolution} * @public (unit-tests) */ solveItWithSubdivisions( dt: number ) { return this.solveWithSubdivisions2( dt ).getFinalState().ltaSolution; } /** * @param {number} dt * @returns {LTACircuit} * @public */ update( dt: number ) { return this.updateCircuit( this.solvePropagate( dt ) ); } /** * Applies the specified solution to the circuit. * * @param {LTASolution} solution * @returns {LTACircuit} * @public */ updateCircuit( solution: LTASolution ) { const updatedCapacitors = this.ltaCapacitors.map( ltaCapacitor => { return new LTACapacitor( ltaCapacitor.id, ltaCapacitor.node0, ltaCapacitor.node1, solution.getVoltage( ltaCapacitor.node0, ltaCapacitor.capacitorVoltageNode1! ), solution.getCurrentForCompanion( ltaCapacitor ), ltaCapacitor.capacitance ); } ); const updatedInductors = this.ltaInductors.map( ltaInductor => { return new LTAInductor( ltaInductor.id, ltaInductor.node0, ltaInductor.node1, solution.getVoltage( ltaInductor.node0, ltaInductor.inductorVoltageNode1! ), solution.getCurrentForCompanion( ltaInductor ), ltaInductor.inductance ); } ); return new LTACircuit( this.ltaResistors, this.ltaBatteries, updatedCapacitors, updatedInductors ); } } /** * @param {number[]} x * @param {number[]} y * @returns {number} */ const euclideanDistance = ( x: number[], y: number[] ) => { assert && assert( x.length === y.length, 'Vector length mismatch' ); let sumSqDiffs = 0; for ( let i = 0; i < x.length; i++ ) { sumSqDiffs += Math.pow( x[ i ] - y[ i ], 2 ); } return Math.sqrt( sumSqDiffs ); }; circuitConstructionKitCommon.register( 'LTACircuit', LTACircuit ); export default LTACircuit; ```

And it is giving pretty good results in the unit test.

Here is the full patch with the forward euler norton inductor:

```diff Index: js/model/analysis/LTAInductor.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/LTAInductor.ts b/js/model/analysis/LTAInductor.ts --- a/js/model/analysis/LTAInductor.ts (revision 7d0c235bc571a01ee979f8cab43cd31a6f05e952) +++ b/js/model/analysis/LTAInductor.ts (date 1637118445093) @@ -7,6 +7,8 @@ constructor( id: number, node0: string, node1: string, voltage: number, current: number, inductance: number ) { super( id, node0, node1, voltage, current ); + + assert && assert( !isNaN( inductance ), 'inductance cannot be NaN' ); this.inductance = inductance; // Synthetic node to read the voltage different across the inductor part (since it is modeled in series with a resistor) Index: js/model/analysis/LTACircuitTests.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/LTACircuitTests.ts b/js/model/analysis/LTACircuitTests.ts --- a/js/model/analysis/LTACircuitTests.ts (revision 7d0c235bc571a01ee979f8cab43cd31a6f05e952) +++ b/js/model/analysis/LTACircuitTests.ts (date 1637106059763) @@ -32,10 +32,11 @@ const t = i * dt; const companionSolution = circuit.solveItWithSubdivisions( dt ); - const voltage = companionSolution!.getVoltage( resistor.nodeId0, resistor.nodeId1 ); - const desiredVoltageAtTPlusDT = -v * Math.exp( -( t + dt ) / r / c ); - const error = Math.abs( voltage - desiredVoltageAtTPlusDT ); + const actualVoltage = companionSolution!.getVoltage( resistor.nodeId0, resistor.nodeId1 ); + const expectedVoltage = v * Math.exp( -( t + dt ) / r / c ); + const error = Math.abs( actualVoltage - expectedVoltage ); + // console.log( expectedVoltage, actualVoltage ); // console.log( error ); if ( error > worstError ) { worstError = error; @@ -45,7 +46,13 @@ // @ts-ignore // window.string = window.string + `${v} ${r} ${c} ${t} ${desiredVoltageAtTPlusDT} ${voltage} ${error} // `; - assert.ok( error < errorThreshold ); // sample run indicates largest error is 1.5328E-7 + + const fractionalError = error / expectedVoltage; + // console.log( fractionalError, error ); + if ( error > 1E-8 ) { + assert.ok( fractionalError <= 0.02 ); // sample run indicates largest error is 1.5328E-7 + } + circuit = circuit.updateWithSubdivisions( dt ); } }; @@ -127,10 +134,13 @@ for ( let i = 0; i < ITERATIONS; i++ ) { const t = i * dt; const solution = circuit.solveItWithSubdivisions( dt ); - const current = solution!.getCurrent( resistor )!; - const expectedCurrent = V / R * ( 1 - Math.exp( -( t + dt ) * R / L ) );//positive, by definition of MNA.Battery - const error = Math.abs( current - expectedCurrent ); - assert.ok( error < errorThreshold ); + const actualCurrent = solution!.getCurrent( resistor )!; + const expectedCurrent = -V / R * ( 1 - Math.exp( -( t + dt ) * R / L ) );//positive, by definition of MNA.Battery + + console.log( expectedCurrent, actualCurrent ); + const error = Math.abs( actualCurrent - expectedCurrent ); + // assert.ok( error < errorThreshold ); + assert.ok( true ); circuit = circuit.updateWithSubdivisions( dt ); } }; @@ -146,8 +156,8 @@ const Leq = L1 + L2; const resistor = new MNAResistor( '1', '2', R ); const battery = new LTAResistiveBattery( id++, '0', '1', V, 0 ); - const inductor1 = new LTAInductor( id++, '2', '3', V, 0.0, L1 ); - const inductor2 = new LTAInductor( id++, '3', '0', V, 0.0, L2 ); + const inductor1 = new LTAInductor( id++, '2', '3', 0, 0.0, L1 ); + const inductor2 = new LTAInductor( id++, '3', '0', 0, 0.0, L2 ); const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor1, inductor2 ] ); iterateInductor( circuit, resistor, V, R, Leq, assert ); }; @@ -155,29 +165,29 @@ const Leq = 1 / ( 1 / L1 + 1 / L2 ); const resistor = new MNAResistor( '1', '2', R ); const battery = new LTAResistiveBattery( id++, '0', '1', V, 0 ); - const inductor1 = new LTAInductor( id++, '2', '0', V, 0.0, L1 ); - const inductor2 = new LTAInductor( id++, '2', '0', V, 0.0, L2 ); + const inductor1 = new LTAInductor( id++, '2', '0', 0, 0.0, L1 ); + const inductor2 = new LTAInductor( id++, '2', '0', 0, 0.0, L2 ); const circuit = new LTACircuit( [ resistor ], [ battery ], [], [ inductor1, inductor2 ] ); iterateInductor( circuit, resistor, V, R, Leq, assert ); }; QUnit.test( 'test_RL_Circuit_should_have_correct_behavior', assert => { testVRLCircuit( 5, 10, 1, assert ); - testVRLCircuit( 3, 11, 2.5, assert ); - testVRLCircuit( 7, 13, 1E4, assert ); - testVRLCircuit( 7, 13, 1E-1, assert ); + // testVRLCircuit( 3, 11, 2.5, assert ); + // testVRLCircuit( 7, 13, 1E4, assert ); + // testVRLCircuit( 7, 13, 1E-1, assert ); } ); - -QUnit.test( 'Series inductors', assert => { - testVRLCircuitSeries( 5, 10, 1, 1, assert ); - for ( let i = 0; i < 10; i++ ) { - testVRLCircuitSeries( 10, 10, 5 * Math.random() + 0.1, 5 * Math.random() + 0.1, assert )// eslint-disable-line - } -} ); - -QUnit.test( 'Parallel inductors', assert => { - testVRLCircuitParallel( 5, 10, 1, 1, assert ); - for ( let i = 0; i < 10; i++ ) { - testVRLCircuitParallel( 10, 10, 5 * Math.random() + 1, 5 * Math.random() + 1, assert )// eslint-disable-line - } -} ); \ No newline at end of file +// +// QUnit.test( 'Series inductors', assert => { +// testVRLCircuitSeries( 5, 10, 1, 1, assert ); +// for ( let i = 0; i < 10; i++ ) { +// testVRLCircuitSeries( 10, 10, 5 * Math.random() + 0.1, 5 * Math.random() + 0.1, assert )// eslint-disable-line +// } +// } ); +// +// QUnit.test( 'Parallel inductors', assert => { +// testVRLCircuitParallel( 5, 10, 1, 1, assert ); +// for ( let i = 0; i < 10; i++ ) { +// testVRLCircuitParallel( 10, 10, 5 * Math.random() + 1, 5 * Math.random() + 1, assert )// eslint-disable-line +// } +// } ); \ No newline at end of file Index: js/model/analysis/mna/MNACircuitTests.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/mna/MNACircuitTests.ts b/js/model/analysis/mna/MNACircuitTests.ts --- a/js/model/analysis/mna/MNACircuitTests.ts (revision 7d0c235bc571a01ee979f8cab43cd31a6f05e952) +++ b/js/model/analysis/mna/MNACircuitTests.ts (date 1637106197461) @@ -12,266 +12,356 @@ import MNABattery from './MNABattery.js'; import MNACurrent from './MNACurrent.js'; import MNACircuitElement from './MNACircuitElement.js'; +import LTAResistiveBattery from '../LTAResistiveBattery.js'; +import LTACircuit from '../LTACircuit.js'; QUnit.module( 'MNACircuitTests' ); const approxEquals = ( a: number, b: number ) => Math.abs( a - b ) < 1E-6; +const approxEquals2 = ( a: number, b: number ) => Math.abs( a - b ) < 1E-4; + +// QUnit.test( 'test_should_be_able_to_obtain_current_for_a_resistor', assert => { +// const battery = new MNABattery( '0', '1', 4.0 ); +// const resistor = new MNAResistor( '1', '0', 2.0 ); +// const solution = new MNACircuit( [ battery ], [ resistor ], [] ).solve(); +// const desiredSolution = new MNASolution( new Map( [ +// [ '0', 0 ], +// [ '1', 4 ] +// ] ), new Map( [ +// [ battery, 2.0 ] +// ] ) ); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solution should match' ); +// +// // same magnitude as battery: positive because current flows from node 1 to 0 +// assert.equal( +// approxEquals( solution.getCurrentForResistor( resistor ), 2 ), true, 'current through resistor should be 2.0 Amps' +// ); +// } ); +// +// QUnit.test( 'test_an_unconnected_resistor_shouldnt_cause_problems', assert => { +// const battery = new MNABattery( '0', '1', 4.0 ); +// const resistor1 = new MNAResistor( '1', '0', 4.0 ); +// const resistor2 = new MNAResistor( '2', '3', 100 ); +// const circuit = new MNACircuit( [ battery ], [ resistor1, resistor2 ], [] ); +// const desiredSolution = new MNASolution( new Map( [ +// [ '0', 0 ], +// [ '1', 4 ], +// [ '2', 0 ], +// [ '3', 0 ] +// ] ), new Map( [ +// [ battery, 1.0 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// -QUnit.test( 'test_battery_resistor_circuit_should_have_correct_voltages_and_currents_for_a_simple_circuit', - assert => { - const battery = new MNABattery( '0', '1', 4.0 ); - const resistor = new MNAResistor( '1', '0', 4.0 ); - const circuit = new MNACircuit( [ battery ], [ resistor ], [] ); +// +// QUnit.test( 'test_current_should_be_reversed_when_voltage_is_reversed', assert => { +// const battery = new MNABattery( '0', '1', -4 ); +// const resistor = new MNAResistor( '1', '0', 2 ); +// const circuit = new MNACircuit( [ battery ], [ resistor ], [] ); +// const voltageMap = new Map( [ +// [ '0', 0 ], +// [ '1', -4 ] +// ] ); +// +// const desiredSolution = new MNASolution( voltageMap, new Map( [ +// [ battery, -2.0 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// +// QUnit.test( 'test_two_batteries_in_series_should_have_voltage_added', assert => { +// const battery1 = new MNABattery( '0', '1', -4 ); +// const battery2 = new MNABattery( '1', '2', -4 ); +// const resistor1 = new MNAResistor( '2', '0', 2.0 ); +// const circuit = new MNACircuit( [ battery1, battery2 ], [ resistor1 ], [] ); +// +// const desiredSolution = new MNASolution( new Map( [ +// [ '0', 0 ], +// [ '1', -4 ], +// [ '2', -8 ] +// ] ), new Map( [ +// [ battery1, -4.0 ], +// [ battery2, -4.0 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// +// QUnit.test( 'test_two_resistors_in_series_should_have_resistance_added', assert => { +// const battery = new MNABattery( '0', '1', 5.0 ); +// const resistor1 = new MNAResistor( '1', '2', 10.0 ); +// const resistor2 = new MNAResistor( '2', '0', 10.0 ); +// const circuit = new MNACircuit( [ battery ], [ +// resistor1, +// resistor2 +// ], [] ); +// const voltageMap = new Map( [ +// [ '0', 0 ], +// [ '1', 5 ], +// [ '2', 2.5 ] +// ] ); +// const desiredSolution = new MNASolution( voltageMap, new Map( [ +// [ battery, 5 / 20.0 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// +// QUnit.test( 'test_A_resistor_with_one_node_unconnected_shouldnt_cause_problems', assert => { +// const battery = new MNABattery( '0', '1', 4.0 ); +// const resistor1 = new MNAResistor( '1', '0', 4.0 ); +// const resistor2 = new MNAResistor( '0', '2', 100.0 ); +// const circuit = new MNACircuit( +// [ battery ], +// [ resistor1, resistor2 ], [] +// ); +// const voltageMap = new Map( [ +// [ '0', 0 ], +// [ '1', 4 ], +// [ '2', 0 ] +// ] ); +// const desiredSolution = new MNASolution( voltageMap, new Map( [ +// [ battery, 1.0 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// +// QUnit.test( 'test_an_unconnected_resistor_shouldnt_cause_problems', assert => { +// const battery = new MNABattery( '0', '1', 4.0 ); +// const resistor1 = new MNAResistor( '1', '0', 4.0 ); +// const resistor2 = new MNAResistor( '2', '3', 100.0 ); +// const circuit = new MNACircuit( [ battery ], [ +// resistor1, +// resistor2 +// ], [] ); +// const voltageMap = new Map( [ +// [ '0', 0 ], +// [ '1', 4 ], +// [ '2', 0 ], +// [ '3', 0 ] +// ] ); +// +// const desiredSolution = new MNASolution( voltageMap, new Map( [ +// [ battery, 1.0 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// +// QUnit.test( 'test_should_handle_resistors_with_no_resistance', assert => { +// const battery = new MNABattery( '0', '1', 5 ); +// const resistor = new MNAResistor( '2', '0', 0 ); +// const resistor0 = new MNAResistor( '1', '2', 10 ); +// const circuit = new MNACircuit( [ battery ], [ +// resistor0, +// resistor +// ], [] ); +// const voltageMap = new Map( [ +// [ '0', 0 ], +// [ '1', 5 ], +// [ '2', 0 ] +// ] ); +// const desiredSolution = new MNASolution( voltageMap, new Map( [ +// [ battery, 5 / 10 ], +// [ resistor, 5 / 10 ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// +// QUnit.test( 'test_resistors_in_parallel_should_have_harmonic_mean_of_resistance', assert => { +// const V = 9.0; +// const R1 = 5.0; +// const R2 = 5.0; +// const Req = 1 / ( 1 / R1 + 1 / R2 ); +// const battery = new MNABattery( '0', '1', V ); +// const resistor1 = new MNAResistor( '1', '0', R1 ); +// const resistor2 = new MNAResistor( '1', '0', R2 ); +// const circuit = new MNACircuit( [ battery ], [ +// resistor1, +// resistor2 +// ], [] ); +// const voltageMap = new Map( [ +// [ '0', 0 ], +// [ '1', V ] +// ] ); +// +// const desiredSolution = new MNASolution( voltageMap, new Map( [ +// [ battery, V / Req ] +// ] ) ); +// const solution = circuit.solve(); +// assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); +// } ); +// - const desiredSolution = new MNASolution( new Map( [ [ '0', 0 ], [ '1', 4 ] ] ), new Map( [ - [ battery, 1.0 ] - ] ) ); - const solution = circuit.solve(); - assert.equal( true, solution.approxEquals( desiredSolution, assert ), 'solutions instances should match' ); - const currentThroughResistor = solution.getCurrentForResistor( resistor ); +// ( () => { +// +// // http://zacg.github.io/blog/2013/08/02/binary-combinations-in-javascript/ +// function binaryCombos( n: number ) { +// const result = []; +// for ( let y = 0; y < Math.pow( 2, n ); y++ ) { +// const combo = []; +// for ( let x = 0; x < n; x++ ) { +// //shift bit and and it with 1 +// if ( ( y >> x ) & 1 ) { // eslint-disable-line +// combo.push( true ); +// } +// else { +// combo.push( false ); +// } +// } +// result.push( combo ); +// } +// return result; +// } +// +// //Usage +// const combos = binaryCombos( 7 ); +// console.log( 'testing ' + combos.length + ' combos' ); +// +// for ( let x = 0; x < combos.length; x++ ) { +// // console.log( combos[ x ].join( ',' ) ); +// +// // @ts-ignore +// window.signs = combos[ x ].map( t => t ? 1 : -1 ); +// +// const v1 = new MNABattery( '1', '0', 24 ); +// const v2 = new MNABattery( '3', '0', 15 ); +// const r1 = new MNAResistor( '1', '2', 10000 ); +// const r2 = new MNAResistor( '2', '3', 8100 ); +// const r3 = new MNAResistor( '2', '0', 4700 ); +// +// const circuit = new MNACircuit( [ v1, v2 ], [ r1, r2, r3 ], [] ); +// const solution = circuit.solve(); +// +// // console.log( solution.getNodeVoltage( '0' ) ); +// // console.log( solution.getNodeVoltage( '1' ) ); +// // console.log( solution.getNodeVoltage( '2' ) ); +// // console.log( solution.getNodeVoltage( '3' ) ); +// // console.log( solution.getSolvedCurrent( v1 ) ); +// // console.log( solution.getSolvedCurrent( v2 ) ); +// +// let wins = 0; +// if ( approxEquals2( solution.getNodeVoltage( '0' ), 0 ) ) { wins++; } +// if ( approxEquals2( solution.getNodeVoltage( '1' ), 24 ) ) { wins++; } +// if ( approxEquals2( solution.getNodeVoltage( '2' ), 9.7470 ) ) { wins++; } +// if ( approxEquals2( solution.getNodeVoltage( '3' ), 15.0 ) ) { wins++; } +// if ( approxEquals2( solution.getSolvedCurrent( v1 ), -1.425E-3 ) ) { wins++; } +// if ( approxEquals2( solution.getSolvedCurrent( v2 ), -6.485E-4 ) ) { wins++; } +// +// +// console.log( 'wins: ' + wins ); +// if ( wins === 6 ) { +// // @ts-ignore +// console.log( window.signs ); +// } +// } +// +// // v1 1 0 dc 24V +// // v2 3 0 dc 15V +// // r1 1 2 10000 Ohms +// // r2 2 3 8100 Ohms +// // r3 2 0 4700 Ohms +// +// // node voltages: (1) 24.0000 (2) 9.7470 (3) 15.0000 +// // voltage source currents: (v1) -1.425E-03 (v2) -6.485E-04 +// +// } )(); - // should be flowing forward through resistor - assert.equal( approxEquals( currentThroughResistor, 1.0 ), true, 'current should be 1 amp through the resistor' ); - } ); +// Netlist example at https://www.allaboutcircuits.com/textbook/reference/chpt-7/example-circuits-and-netlists/ +QUnit.test( 'Netlist: 2 batteries & 3 resistors', assert => { -QUnit.test( 'test_battery_resistor_circuit_should_have_correct_voltages_and_currents_for_a_simple_circuit_ii', - assert => { - const battery = new MNABattery( '0', '1', 4.0 ); - const resistor = new MNAResistor( '1', '0', 2.0 ); - const circuit = new MNACircuit( [ battery ], [ resistor ], [] ); - const desiredSolution = new MNASolution( new Map( [ - [ '0', 0 ], - [ '1', 4 ] - ] ), new Map( [ - [ battery, 2.0 ] - ] ) ); - const solution = circuit.solve(); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solution should match' ); - } ); - - -QUnit.test( 'test_should_be_able_to_obtain_current_for_a_resistor', assert => { - const battery = new MNABattery( '0', '1', 4.0 ); - const resistor = new MNAResistor( '1', '0', 2.0 ); - const solution = new MNACircuit( [ battery ], [ resistor ], [] ).solve(); - const desiredSolution = new MNASolution( new Map( [ - [ '0', 0 ], - [ '1', 4 ] - ] ), new Map( [ - [ battery, 2.0 ] + const v1 = new MNABattery( '1', '0', 24 ); + const v2 = new MNABattery( '3', '0', 15 ); + const r1 = new MNAResistor( '1', '2', 10000 ); + const r2 = new MNAResistor( '2', '3', 8100 ); + const r3 = new MNAResistor( '2', '0', 4700 ); + + const circuit = new MNACircuit( [ v1, v2 ], [ r1, r2, r3 ], [] ); + + const desiredSolution = new MNASolution( new Map( [ + [ '0', 0 ], + [ '1', 24.0 ], + [ '2', 9.7470 ], + [ '3', 15.0 ] + ] ), new Map( [ + [ v1, -1.425E-3 ], + [ v2, -6.485E-4 ] ] ) ); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solution should match' ); - - // same magnitude as battery: positive because current flows from node 1 to 0 - assert.equal( - approxEquals( solution.getCurrentForResistor( resistor ), 2 ), true, 'current through resistor should be 2.0 Amps' - ); -} ); - -QUnit.test( 'test_an_unconnected_resistor_shouldnt_cause_problems', assert => { - const battery = new MNABattery( '0', '1', 4.0 ); - const resistor1 = new MNAResistor( '1', '0', 4.0 ); - const resistor2 = new MNAResistor( '2', '3', 100 ); - const circuit = new MNACircuit( [ battery ], [ resistor1, resistor2 ], [] ); - const desiredSolution = new MNASolution( new Map( [ - [ '0', 0 ], - [ '1', 4 ], - [ '2', 0 ], - [ '3', 0 ] - ] ), new Map( [ - [ battery, 1.0 ] - ] ) ); - const solution = circuit.solve(); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); -} ); - -QUnit.test( 'test_current_source_should_provide_current', assert => { - const currentSource = new MNACurrent( '0', '1', 10 ); - const resistor = new MNAResistor( '1', '0', 4 ); - const circuit = new MNACircuit( [], [ resistor ], [ currentSource ] ); - const voltageMap = new Map( [ - [ '0', 0 ], - - // This is negative since traversing across the resistor should yield a negative voltage, see - // http://en.wikipedia.org/wiki/Current_source - [ '1', -40.0 ] - ] ); - const desiredSolution = new MNASolution( voltageMap, new Map() ); const solution = circuit.solve(); + assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); } ); -QUnit.test( 'test_current_should_be_reversed_when_voltage_is_reversed', assert => { - const battery = new MNABattery( '0', '1', -4 ); - const resistor = new MNAResistor( '1', '0', 2 ); - const circuit = new MNACircuit( [ battery ], [ resistor ], [] ); - const voltageMap = new Map( [ - [ '0', 0 ], - [ '1', -4 ] - ] ); - - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, -2.0 ] - ] ) ); +QUnit.test( 'Netlist: 1 battery and 1 resistor', assert => { + const v1 = new MNABattery( 'a', '0', 9 ); + const r1 = new MNAResistor( 'a', '0', 9 ); + const circuit = new MNACircuit( [ v1 ], [ r1 ], [] ); + const desiredSolution = new MNASolution( new Map( [ [ 'a', 9 ], [ '0', 0 ] ] ), new Map( [ [ v1, -1.0 ] ] ) ); const solution = circuit.solve(); assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); } ); -QUnit.test( 'test_two_batteries_in_series_should_have_voltage_added', assert => { - const battery1 = new MNABattery( '0', '1', -4 ); - const battery2 = new MNABattery( '1', '2', -4 ); - const resistor1 = new MNAResistor( '2', '0', 2.0 ); - const circuit = new MNACircuit( [ battery1, battery2 ], [ resistor1 ], [] ); - - const desiredSolution = new MNASolution( new Map( [ - [ '0', 0 ], - [ '1', -4 ], - [ '2', -8 ] - ] ), new Map( [ - [ battery1, -4.0 ], - [ battery2, -4.0 ] - ] ) ); +QUnit.test( 'Netlist: 1 battery and 1 resistor: part 2', assert => { + const v1 = new MNABattery( 'a', '0', 10 ); + const r1 = new MNAResistor( 'a', '0', 5 ); + const circuit = new MNACircuit( [ v1 ], [ r1 ], [] ); + const desiredSolution = new MNASolution( new Map( [ [ 'a', 10 ], [ '0', 0 ] ] ), new Map( [ [ v1, -2.0 ] ] ) ); const solution = circuit.solve(); assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); } ); -QUnit.test( 'test_two_resistors_in_series_should_have_resistance_added', assert => { - const battery = new MNABattery( '0', '1', 5.0 ); - const resistor1 = new MNAResistor( '1', '2', 10.0 ); - const resistor2 = new MNAResistor( '2', '0', 10.0 ); - const circuit = new MNACircuit( [ battery ], [ - resistor1, - resistor2 - ], [] ); +// Compare to the example at https://www.khanacademy.org/science/electrical-engineering/ee-circuit-analysis-topic/ee-dc-circuit-analysis/a/ee-node-voltage-method +QUnit.test( 'Netlist: khan ee', assert => { + const battery = new MNABattery( 'a', 'c', 140 ); + const resistor1 = new MNAResistor( 'a', 'b', 20 ); + const resistor2 = new MNAResistor( 'b', 'c', 6 ); + const resistor3 = new MNAResistor( 'b', 'c', 5 ); + const current = new MNACurrent( 'c', 'b', 18 ); + const circuit = new MNACircuit( [ battery ], [ resistor1, resistor2, resistor3 ], [ current ] ); + + const offset = -140; const voltageMap = new Map( [ - [ '0', 0 ], - [ '1', 5 ], - [ '2', 2.5 ] + [ 'a', 140 + offset ], + [ 'b', 60 + offset ], + [ 'c', 0 + offset ] ] ); - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, 5 / 20.0 ] + + const desiredSolution = new MNASolution( voltageMap, new Map( [ + [ battery, -4.0 ] ] ) ); const solution = circuit.solve(); assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); } ); -QUnit.test( 'test_A_resistor_with_one_node_unconnected_shouldnt_cause_problems', assert => { - const battery = new MNABattery( '0', '1', 4.0 ); - const resistor1 = new MNAResistor( '1', '0', 4.0 ); - const resistor2 = new MNAResistor( '0', '2', 100.0 ); - const circuit = new MNACircuit( - [ battery ], - [ resistor1, resistor2 ], [] - ); +QUnit.test( 'netlist: ir', assert => { + const currentSource = new MNACurrent( '0', 'a', 10 ); + const resistor = new MNAResistor( 'a', '0', 4 ); + const circuit = new MNACircuit( [], [ resistor ], [ currentSource ] ); const voltageMap = new Map( [ [ '0', 0 ], - [ '1', 4 ], - [ '2', 0 ] + [ 'a', 40.0 ] ] ); - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, 1.0 ] - ] ) ); + const desiredSolution = new MNASolution( voltageMap, new Map() ); const solution = circuit.solve(); assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); } ); -QUnit.test( 'test_an_unconnected_resistor_shouldnt_cause_problems', assert => { - const battery = new MNABattery( '0', '1', 4.0 ); - const resistor1 = new MNAResistor( '1', '0', 4.0 ); - const resistor2 = new MNAResistor( '2', '3', 100.0 ); - const circuit = new MNACircuit( [ battery ], [ - resistor1, - resistor2 - ], [] ); - const voltageMap = new Map( [ - [ '0', 0 ], - [ '1', 4 ], - [ '2', 0 ], - [ '3', 0 ] - ] ); - - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, 1.0 ] - ] ) ); - const solution = circuit.solve(); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); -} ); - -QUnit.test( 'test_should_handle_resistors_with_no_resistance', assert => { - const battery = new MNABattery( '0', '1', 5 ); - const resistor = new MNAResistor( '2', '0', 0 ); - const resistor0 = new MNAResistor( '1', '2', 10 ); - const circuit = new MNACircuit( [ battery ], [ - resistor0, - resistor - ], [] ); - const voltageMap = new Map( [ - [ '0', 0 ], - [ '1', 5 ], - [ '2', 0 ] - ] ); - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, 5 / 10 ], - [ resistor, 5 / 10 ] - ] ) ); - const solution = circuit.solve(); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); -} ); - -QUnit.test( 'test_resistors_in_parallel_should_have_harmonic_mean_of_resistance', assert => { - const V = 9.0; - const R1 = 5.0; - const R2 = 5.0; - const Req = 1 / ( 1 / R1 + 1 / R2 ); - const battery = new MNABattery( '0', '1', V ); - const resistor1 = new MNAResistor( '1', '0', R1 ); - const resistor2 = new MNAResistor( '1', '0', R2 ); - const circuit = new MNACircuit( [ battery ], [ - resistor1, - resistor2 - ], [] ); - const voltageMap = new Map( [ - [ '0', 0 ], - [ '1', V ] - ] ); - - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, V / Req ] - ] ) ); - const solution = circuit.solve(); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); -} ); - -// Compare to the example at https://www.khanacademy.org/science/electrical-engineering/ee-circuit-analysis-topic/ee-dc-circuit-analysis/a/ee-node-voltage-method -QUnit.test( 'test-circuit-with-voltage-and-current-sources', assert => { - const V1 = 140; - const R1 = 20; - const R2 = 6; - const R3 = 5; - const I1 = 18; - - // TODO: Sign error - const battery = new MNABattery( 'a', 'c', -V1 ); - const resistor1 = new MNAResistor( 'a', 'b', R1 ); - const resistor2 = new MNAResistor( 'b', 'c', R2 ); - const resistor3 = new MNAResistor( 'b', 'c', R3 ); - - // TODO: Sign error - const current = new MNACurrent( 'c', 'b', -I1 ); - const circuit = new MNACircuit( [ battery ], [ resistor1, resistor2, resistor3 ], [ current ] ); - - const voltageMap = new Map( [ - [ 'a', 140 - 140 ], - [ 'b', 60 - 140 ], - [ 'c', 0 - 140 ] - ] ); - - const desiredSolution = new MNASolution( voltageMap, new Map( [ - [ battery, 4.0 ] - ] ) ); - const solution = circuit.solve(); - - console.log( solution ); - assert.equal( solution.approxEquals( desiredSolution, assert ), true, 'solutions should match' ); -} ); +// const findOperatingPointForLRCircuit = ( V: number, R: number, L: number ) => { +// let id = 0; +// const resistor = new MNAResistor( '1', '2', R ); +// const battery = new MNABattery( '0', '1', V ); +// const myResistor = new MNAResistor( '2', '0', 1 ); +// const circuit = new MNACircuit( [ battery ], [ resistor, myResistor ], [] ); +// // iterateInductor( circuit, resistor, V, R, L, assert ); +// const x = circuit.solve(); +// debugger; +// }; +// +// findOperatingPointForLRCircuit( 5, 10, 1 ); \ No newline at end of file Index: js/model/analysis/LTACircuit.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/analysis/LTACircuit.ts b/js/model/analysis/LTACircuit.ts --- a/js/model/analysis/LTACircuit.ts (revision 7d0c235bc571a01ee979f8cab43cd31a6f05e952) +++ b/js/model/analysis/LTACircuit.ts (date 1637117551950) @@ -55,6 +55,7 @@ const companionBatteries: MNABattery[] = []; const companionResistors: MNAResistor[] = []; + const companionCurrents: MNACurrent[] = []; const currentCompanions: { element: CoreModel, getValueForSolution: ( solution: MNASolution ) => number }[] = []; // Node indices that have been used @@ -118,40 +119,67 @@ } ); } ); + const signs = [ 1, 1, 1, 1 ]; + // const signs = [ -1, +1, +1, +1 ]; + // const signs = [ +1, -1, +1, +1 ]; + // const signs = [ +1, 1, -1, +1 ]; + // const signs = [ +1, 1, 1, -1 ]; + // const signs = [ -1, -1, 1, 1 ]; + // const signs = [ -1, 1, -1, 1 ]; + // const signs = [ -1, 1, 1, -1 ]; + // const signs = [ 1, -1, -1, 1 ]; + // const signs = [ 1, -1, 1, -1 ]; + // const signs = [ 1, 1, -1, -1 ]; + // const signs = [ 1, -1, -1, -1 ]; + // const signs = [ -1, 1, -1, -1 ]; + // const signs = [ -1, -1, 1, -1 ]; + // const signs = [ -1, -1, -1, 1 ]; + // const signs = [ -1, -1, -1, -1 ]; + // const A = +1; + // const B = +1; + // const C = +1; + // const D = +1; + // See also http://circsimproj.blogspot.com/2009/07/companion-models.html, which reports: // Req = 2L/dt // Veq = -2Li/dt-v // See najm page 279 and Pillage page 86 this.ltaInductors.forEach( ltaInductor => { - // In series - const newNode = 'syntheticNode' + ( syntheticNodeIndex++ ); - const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); + // debugger; - const companionResistance = 2 * ltaInductor.inductance / dt; - const companionVoltage = ltaInductor.voltage + -companionResistance * ltaInductor.current; + // const newNode1 = 'syntheticNode' + ( syntheticNodeIndex++ ); + // const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ ); - const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); - const resistor = new MNAResistor( newNode, newNode2, companionResistance ); - const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); - companionBatteries.push( battery ); - companionResistors.push( resistor ); - companionResistors.push( addedResistor ); + // Req = dt/2/L + // Ieq = i0 + Req * v0 + // const companionResistance = dt / ltaInductor.inductance; + const companionCurrent = -dt * ltaInductor.voltage / ltaInductor.inductance + ltaInductor.current; // signs[ 0 ] * ltaInductor.current + signs[ 1 ] * companionResistance * ltaInductor.voltage; + // console.log( companionCurrent ); - ltaInductor.inductorVoltageNode1 = newNode2; + // const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage ); + // const resistor = new MNAResistor( ltaInductor.node0, ltaInductor.node1, companionResistance ); + companionCurrents.push( new MNACurrent( ltaInductor.node0, ltaInductor.node1, companionCurrent ) ); + // const addedResistor = new MNAResistor( newNode2, ltaInductor.node1, CCKCQueryParameters.inductorResistance ); + // companionBatteries.push( battery ); + // companionResistors.push( resistor ); - // we need to be able to get the current for this component - // in series, so current is same through both companion components + ltaInductor.inductorVoltageNode1 = ltaInductor.node1; + currentCompanions.push( { element: ltaInductor, - getValueForSolution: ( solution: MNASolution ) => solution.getCurrentForResistor( resistor ) + getValueForSolution: ( solution: MNASolution ) => { + return companionCurrent; + + // const a = signs[ 2 ] * solution.getCurrentForResistor( resistor ); + // return a + signs[ 3 ] * companionCurrent; + } } ); } ); const newBatteryList = companionBatteries; const newResistorList = [ ...this.ltaResistors, ...companionResistors ]; - const newCurrentList: MNACurrent[] = []; // Placeholder for if we add other circuit elements in the future - + const newCurrentList = companionCurrents; const mnaCircuit = new MNACircuit( newBatteryList, newResistorList, newCurrentList ); const mnaSolution = mnaCircuit.solve(); ```

The values are coming out like this:

Expected Current	Actual Current
-0.07675913755469299	-0.07852840330803694
-0.14173434471310536	-0.1447233863638541
-0.1967346701436833	-0.20052199676691962
-0.243291440483704	-0.24755705560646554
-0.28270089574646085	-0.2872049383056726
-0.31606027941427883	-0.32062585115905434
-0.3442983880427012	-0.3487977821654892
-0.36820143094213664	-0.37211561741076454
-0.3884349199257851	-0.39183746412946363
-0.40556219858121906	-0.4088251266087866
-0.42006012696015305	-0.42288577139755573
-0.43233235838169365	-0.4347780366260082
-0.44272057800365616	-0.44483632160429243
-0.45151401606779745	-0.45334345584330443
-0.4589575006880506	-0.460538651954452
-0.46525827438859924	-0.46662423208324116
-0.47059176417878507	-0.4717713170176646
-0.475106465816068	-0.4761246379437748
-0.4789280782453618	-0.4798160412860006
-0.4821630033263738	-0.4829367115602657
-0.48490130828884076	-0.48557489060975845
-0.4872192333967463	-0.48780517708205207
-0.48918131464025344	-0.4896906358227885
-0.4908421805556329	-0.4912845811330362
-0.49224807320049535	-0.49263208431471084
-0.4934381356315295	-0.49377124813228396
-0.49444550173087887	-0.4947394317149829
-0.4952982187242524	-0.4955571229567258
-0.49602002807564677	-0.496247714058219
-0.49663102650045726	-0.4968309611875923
-0.4971482255009963	-0.49732354965736497

I also tried putting an exponential damping term on the inductor voltage, and that seemed to correct some bad behavior:

Index: js/model/analysis/LTACircuit.ts
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/js/model/analysis/LTACircuit.ts b/js/model/analysis/LTACircuit.ts
--- a/js/model/analysis/LTACircuit.ts   (revision 0564d2ce40170627e3585b0438c57fe6664bf38f)
+++ b/js/model/analysis/LTACircuit.ts   (date 1637120966663)
@@ -133,7 +133,7 @@
         const newNode2 = 'syntheticNode' + ( syntheticNodeIndex++ );

         const companionResistance = 2 * ltaInductor.inductance / dt;
-        const companionVoltage = ltaInductor.voltage + -companionResistance * ltaInductor.current;
+        const companionVoltage = ( ltaInductor.voltage + -companionResistance * ltaInductor.current ) * 0.99;

         const battery = new MNABattery( ltaInductor.node0, newNode, companionVoltage );
         const resistor = new MNAResistor( newNode, newNode2, companionResistance );

It fixes the bad oscillation in this case:

Unit tests are all still within tolerance of 2% after this damping. 0.999 damping wasn't enough. 0.99 was enough. I didn't search between those values.

After the commit, the entire forward euler norton inductor is this (does not include the 1E-4 series resistance):

    this.ltaInductors.forEach( ltaInductor => {

      const companionCurrent = -dt * ltaInductor.voltage / ltaInductor.inductance + ltaInductor.current; // signs[ 0 ] * ltaInductor.current + signs[ 1 ] * companionResistance * ltaInductor.voltage;
      companionCurrents.push( new MNACurrent( ltaInductor.node0, ltaInductor.node1, companionCurrent ) );
      ltaInductor.inductorVoltageNode1 = ltaInductor.node1;

      currentCompanions.push( {
        element: ltaInductor,
        getValueForSolution: ( solution: MNASolution ) => companionCurrent
      } );
    } );

Today, we observed that putting dt low enough, say to 1E-5, would correct the problem in the case described in https://github.com/phetsims/circuit-construction-kit-common/issues/780#issuecomment-971149638. We did not get the damping working.

However, we aren't sure how to pick a threshold for various architectures and lower inductance values. We observed that lower inductance values fail earlier. We also saw that inductors in Java go from 10 H to 100 H. In HTML5, we have been targeting between 0.1 H and 10 H. The question is how we can get resonance if we have the higher inductance values.

The resonant frequency is sqrt(1/LC) and will need to match the driving frequency.

We want to hit all our learning goals about resonance, and weren’t sure what was possible. You can still use the same query parameters for changing the inductor min/max/default/step. And we added a decimal place (in master) to the ac frequency control for investigation.

At today's design meeting, we discussed that the range in java was 10H to 100H.

@ariel-phet: If you want to hit resonance at 10H, you need to reduce the capacitance accordingly. The capacitance ranges were already finely tuned. @kathy-phet: We need at least 1 way to hit resonance, but do not need to cover a large state space with resonances. Let's calculate where the resonance is. @ariel-phet: Can you please add capacitance range query parameters like we have for inductors? @samreid: Sure. @ariel-phet: There is a sweet spot for resonant frequency, we need to target that region. @arouinfar: Changing the capacitance will change the rate at which the plates saturate. @kathy-phet: I just want to understand if there are too many constraints to hit all of these different spaces. @ariel-phet: We really need to be able to hit resonance. @arouinfar: Maybe we will leave the default for capacitance as we have it now, but add a lower range for the resonance learning goals. @samreid: That sounds good to me. @ariel-phet: What about a screen that is all about resonance? Discussion about the plate distance. @ariel-phet so the proposal is to let the capacitance go to a lower range. Let's investigate that. @samreid: I'll add the query parameters to adjust those values. @kathy-phet : If L=5, we need C=0.005 to get f=1. Currently the min capacitance is 0.05. So it's just a factor of 10 less. @samreid: We will have to make sure the low capacitance doesn't give integration errors. @kathy-phet: Also, if resonance is the only feature we cannot solve soon (like within a week), then we may need to move forward without it, since so many other things are already working much better than before. @ariel-phet will also check Java to see what resonance was possible.

Next step is for @samreid to add capacitance range query parameters.

I added the query parameters, they will be ready for testing in the next dev version:

  capacitanceMin: {
    defaultValue: 0.05
  },
  capacitanceMax: {
    defaultValue: 0.2
  },
  capacitanceStep: {
    defaultValue: 0.001
  },
  capacitanceDefault: {
    defaultValue: 0.1
  },
  capacitorNumberDecimalPlaces: {
    defaultValue: 2
  },

UPDATE: Please test in https://phet-dev.colorado.edu/html/circuit-construction-kit-ac/1.0.0-dev.51/phet/circuit-construction-kit-ac_en_phet.html

I thought after the last group meeting I attended to discuss this, we decided to stick with the original values @ariel-phet had chosen and then restrict to one inductor. So you can hit resonance with the L values Ariel chose a while back.

I thought any revisiting of values was going to happen - as needed - after publishing 1.0, when we continue to diagnose and figure out if we can support a version with multiple inductors allowed.

So, bottom line, my recommendation is to 1) Keep these query parameters in so that we can investigate 2) But not to make this a blocking issue. Use the values that Ariel had found for L that allow access to resonance for now, and the values of R and C that we have had for a long time now.

@arouinfar - Do you agree with this recollection (see ^^)

@kathy-phet that was my understanding, that since we were going to limit to 1 inductor for the moment, no need to mess around with the current default values of L and C if we can avoid it.

@arouinfar and I cleared the milestone since it seems we don't need to change the capacitor ranges or anything else for RC.1. @arouinfar recommended closing this issue.

phetsims / circuit-construction-kit-common

Model lag for inductor capacitor circuit #780