Skip to content

Commit d0dc95a

Browse files
lukelowrylukelowry
committed
clean/correct ieeest impl per review
1 parent 349a26d commit d0dc95a

2 files changed

Lines changed: 95 additions & 52 deletions

File tree

‎GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp‎

Lines changed: 40 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -124,17 +124,18 @@ namespace GridKit
124124
size_ = static_cast<IdxT>(IeeestInternalVariables<order>::MAXIMUM);
125125
auto size = static_cast<size_t>(size_);
126126

127-
f_.assign(size, ScalarT{0});
128-
y_.assign(size, ScalarT{0});
129-
yp_.assign(size, ScalarT{0});
130-
tag_.assign(size, false);
131-
abs_tol_.assign(size, ScalarT{0});
127+
if (!allocated_)
128+
{
129+
this->allocateVectors(size_);
130+
}
131+
132+
tag_.resize(size);
132133
variable_indices_.resize(size);
133134
residual_indices_.resize(size);
134135

135136
auto signal_size = static_cast<size_t>(IeeestExternalVariables::MAXIMUM);
136-
ws_.assign(signal_size, ScalarT{0});
137-
ws_indices_.assign(signal_size, INVALID_INDEX<IdxT>);
137+
ws_.resize(signal_size);
138+
ws_indices_.resize(signal_size);
138139

139140
for (IdxT j = 0; j < size_; ++j)
140141
{
@@ -145,11 +146,13 @@ namespace GridKit
145146
constexpr auto VSS = IeeestInternalVariables<order>::VSS;
146147
if (signals_.template isAssigned<VSS>())
147148
{
149+
auto* y = y_.getData();
148150
signals_.template getSignalNode<VSS>()->set(
149-
&y_[static_cast<size_t>(VSS)],
151+
&y[static_cast<size_t>(VSS)],
150152
&(this->getVariableIndex(static_cast<IdxT>(VSS))));
151153
}
152154

155+
allocated_ = true;
153156
return 0;
154157
}
155158

@@ -220,8 +223,13 @@ namespace GridKit
220223
const auto VSS = static_cast<size_t>(IeeestInternalVariables<order>::VSS);
221224
const auto U = static_cast<size_t>(IeeestExternalVariables::U);
222225

223-
std::fill(y_.begin(), y_.end(), ZERO<RealT>);
224-
std::fill(yp_.begin(), yp_.end(), ZERO<RealT>);
226+
auto* y = y_.getData();
227+
auto* yp = yp_.getData();
228+
for (IdxT i = 0; i < size_; ++i)
229+
{
230+
y[i] = ZERO<RealT>;
231+
yp[i] = ZERO<RealT>;
232+
}
225233

226234
const ScalarT u = signals_.template readExternalVariable<IeeestExternalVariables::U>();
227235
ws_[U] = u;
@@ -233,17 +241,20 @@ namespace GridKit
233241
{
234242
const auto X1 = static_cast<size_t>(IeeestInternalVariables<order>::X1);
235243

236-
y_[X1] = u;
244+
y[X1] = u;
237245
}
238246

239-
y_[X5] = u;
240-
y_[X6] = u;
241-
y_[X7] = u;
242-
y_[V4] = u;
243-
y_[V5] = u;
244-
y_[V6] = u;
245-
y_[V7] = ZERO<RealT>;
246-
y_[VSS] = Math::clamp(y_[V7], Lsmin_, Lsmax_);
247+
y[X5] = u;
248+
y[X6] = u;
249+
y[X7] = u;
250+
y[V4] = u;
251+
y[V5] = u;
252+
y[V6] = u;
253+
y[V7] = ZERO<RealT>;
254+
y[VSS] = Math::clamp(y[V7], Lsmin_, Lsmax_);
255+
256+
y_.setDataUpdated();
257+
yp_.setDataUpdated();
247258

248259
return 0;
249260
}
@@ -292,7 +303,7 @@ namespace GridKit
292303
template <typename scalar_type, typename index_type, size_t order>
293304
int Ieeest<scalar_type, index_type, order>::setAbsoluteTolerance(RealT rel_tol)
294305
{
295-
std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol);
306+
abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
296307
return 0;
297308
}
298309

@@ -428,16 +439,21 @@ namespace GridKit
428439
{
429440
const auto U = static_cast<size_t>(IeeestExternalVariables::U);
430441

431-
std::fill(ws_.begin(), ws_.end(), ZERO<RealT>);
432-
std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX<IdxT>);
442+
ws_[U] = ZERO<RealT>;
443+
ws_indices_[U] = INVALID_INDEX<IdxT>;
433444

434445
if (signals_.template isAttached<IeeestExternalVariables::U>())
435446
{
436447
ws_[U] = signals_.template readExternalVariable<IeeestExternalVariables::U>();
437448
ws_indices_[U] = signals_.template readExternalVariableIndex<IeeestExternalVariables::U>();
438449
}
439450

440-
evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data());
451+
const auto* y = y_.getData();
452+
const auto* yp = yp_.getData();
453+
auto* f = f_.getData();
454+
evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
455+
456+
f_.setDataUpdated();
441457

442458
return 0;
443459
}
@@ -456,7 +472,7 @@ namespace GridKit
456472
constexpr auto VSS = static_cast<size_t>(IeeestInternalVariables<order>::VSS);
457473

458474
monitor_->set(Variable::vss, [this]
459-
{ return y_[VSS]; });
475+
{ return y_.getData()[VSS]; });
460476
}
461477

462478
} // namespace Stabilizer

‎tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp‎

Lines changed: 55 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -111,11 +111,12 @@ namespace GridKit
111111
success *= (model.verify() == 0);
112112
model.initialize();
113113

114-
success *= vss_node.linked();
115-
success *= (vss_node.getVariableIndex() == static_cast<IdxT>(VSS));
116-
success *= isEqual(model.y()[V7], test.expected_v7, tol_);
117-
success *= isEqual(model.y()[VSS], test.expected_vss, loose_tol);
118-
success *= isEqual(vss_node.read(), test.expected_vss, loose_tol);
114+
success *= vss_node.linked();
115+
success *= (vss_node.getVariableIndex() == static_cast<IdxT>(VSS));
116+
const auto* y = model.y().getData();
117+
success *= isEqual(y[V7], test.expected_v7, tol_);
118+
success *= isEqual(y[VSS], test.expected_vss, loose_tol);
119+
success *= isEqual(vss_node.read(), test.expected_vss, loose_tol);
119120
}
120121

121122
const std::string name = orderedName(__func__, order);
@@ -162,13 +163,15 @@ namespace GridKit
162163
success *= (model.evaluateResidual() == 0);
163164

164165
// The smooth clamp keeps the VSS row only approximately zero.
165-
const auto loose_tol = static_cast<RealT>(1.0e-4);
166-
for (size_t i = 0; i < model.getResidual().size(); ++i)
166+
const auto loose_tol = static_cast<RealT>(1.0e-4);
167+
const auto& residual = model.getResidual();
168+
const auto* f = residual.getData();
169+
for (size_t i = 0; i < residual.getSize(); ++i)
167170
{
168-
if (!isEqual(model.getResidual()[i], static_cast<ScalarT>(0.0), loose_tol))
171+
if (!isEqual(f[i], static_cast<ScalarT>(0.0), loose_tol))
169172
{
170173
std::cout << "Nonzero initial residual at row " << i << ": "
171-
<< std::setprecision(15) << model.getResidual()[i] << "\n";
174+
<< std::setprecision(15) << f[i] << "\n";
172175
success = false;
173176
}
174177
}
@@ -202,13 +205,18 @@ namespace GridKit
202205

203206
const auto y_values = stateValues<order>();
204207
const auto yp_values = derivativeValues<order>();
208+
auto* y = model.y().getData();
209+
auto* yp = model.yp().getData();
205210
for (size_t i = 0; i < y_values.size(); ++i)
206211
{
207-
model.y()[i] = y_values[i];
208-
model.yp()[i] = yp_values[i];
212+
y[i] = y_values[i];
213+
yp[i] = yp_values[i];
209214
}
215+
model.y().setDataUpdated();
216+
model.yp().setDataUpdated();
210217

211218
model.evaluateResidual();
219+
const auto* f = model.getResidual().getData();
212220

213221
// The smooth clamp on the VSS row carries approximation error, so
214222
// that row is compared with a looser tolerance.
@@ -218,11 +226,11 @@ namespace GridKit
218226
for (size_t i = 0; i < expected.size(); ++i)
219227
{
220228
const auto test_tol = (i == VSS) ? loose_tol : tol_;
221-
if (!isEqual(model.getResidual()[i], expected[i], test_tol))
229+
if (!isEqual(f[i], expected[i], test_tol))
222230
{
223231
std::cout << "Incorrect residual for order " << order
224232
<< " row " << i << ": "
225-
<< std::setprecision(15) << model.getResidual()[i]
233+
<< std::setprecision(15) << f[i]
226234
<< " != " << expected[i] << "\n";
227235
success = false;
228236
}
@@ -443,13 +451,18 @@ namespace GridKit
443451

444452
const auto y_values = stateValues<4>();
445453
const auto yp_values = derivativeValues<4>();
454+
auto* y = stabilizer->y().getData();
455+
auto* yp = stabilizer->yp().getData();
446456
for (size_t i = 0; i < y_values.size(); ++i)
447457
{
448-
stabilizer->y()[i] = y_values[i];
449-
stabilizer->yp()[i] = yp_values[i];
458+
y[i] = y_values[i];
459+
yp[i] = yp_values[i];
450460
}
461+
stabilizer->y().setDataUpdated();
462+
stabilizer->yp().setDataUpdated();
451463

452464
stabilizer->evaluateResidual();
465+
const auto* f = stabilizer->getResidual().getData();
453466

454467
// Derived with the a1 and a3 residual terms dropped:
455468
// x4_rhs = (0.5 - 0.1 - 0.6 * 0.3) / 0.08 = 2.75
@@ -463,10 +476,10 @@ namespace GridKit
463476
for (size_t i = 0; i < expected.size(); ++i)
464477
{
465478
const auto test_tol = (i == VSS) ? loose_tol : tol_;
466-
if (!isEqual(stabilizer->getResidual()[i], expected[i], test_tol))
479+
if (!isEqual(f[i], expected[i], test_tol))
467480
{
468481
std::cout << "Incorrect symmetric-notch residual row " << i << ": "
469-
<< std::setprecision(15) << stabilizer->getResidual()[i]
482+
<< std::setprecision(15) << f[i]
470483
<< " != " << expected[i] << "\n";
471484
success = false;
472485
}
@@ -510,44 +523,54 @@ namespace GridKit
510523
model.allocate();
511524
model.initialize();
512525

526+
auto* y = model.y().getData();
527+
auto* yp = model.yp().getData();
513528
for (size_t i = 0; i < model.size(); ++i)
514529
{
515-
model.y()[i].setVariableNumber(i);
530+
y[i].setVariableNumber(i);
516531
}
517532
u_value.setVariableNumber(model.size());
518533
u_value.setValue(0.5);
519534

520535
for (size_t i = 0; i < y_values.size(); ++i)
521536
{
522-
model.y()[i].setValue(y_values[i]);
537+
y[i].setValue(y_values[i]);
523538
}
524539
for (size_t i = 0; i < yp_values.size(); ++i)
525540
{
526-
model.yp()[i].setValue(yp_values[i]);
541+
yp[i].setValue(yp_values[i]);
527542
}
543+
model.y().setDataUpdated();
544+
model.yp().setDataUpdated();
528545

529546
model.evaluateResidual();
530-
std::vector<DepVar> residual_y = model.getResidual();
547+
const auto& residual_y_view = model.getResidual();
548+
std::vector<DepVar> residual_y(residual_y_view.getData(),
549+
residual_y_view.getData() + residual_y_view.getSize());
531550

532551
model.initialize();
533552
for (size_t i = 0; i < model.size(); ++i)
534553
{
535-
model.y()[i] = model.y()[i].getValue();
536-
model.yp()[i].setVariableNumber(i);
554+
y[i] = y[i].getValue();
555+
yp[i].setVariableNumber(i);
537556
}
538557
u_value = 0.5;
539558

540559
for (size_t i = 0; i < y_values.size(); ++i)
541560
{
542-
model.y()[i].setValue(y_values[i]);
561+
y[i].setValue(y_values[i]);
543562
}
544563
for (size_t i = 0; i < yp_values.size(); ++i)
545564
{
546-
model.yp()[i].setValue(yp_values[i]);
565+
yp[i].setValue(yp_values[i]);
547566
}
567+
model.y().setDataUpdated();
568+
model.yp().setDataUpdated();
548569

549570
model.evaluateResidual();
550-
std::vector<DepVar> residual_yp = model.getResidual();
571+
const auto& residual_yp_view = model.getResidual();
572+
std::vector<DepVar> residual_yp(residual_yp_view.getData(),
573+
residual_yp_view.getData() + residual_yp_view.getSize());
551574

552575
dependency_tracking_jacobian.resize(residual_y.size());
553576
for (size_t i = 0; i < residual_y.size(); ++i)
@@ -598,11 +621,15 @@ namespace GridKit
598621
model.allocate();
599622
model.initialize();
600623

624+
auto* y = model.y().getData();
625+
auto* yp = model.yp().getData();
601626
for (size_t i = 0; i < y_values.size(); ++i)
602627
{
603-
model.y()[i] = y_values[i];
604-
model.yp()[i] = yp_values[i];
628+
y[i] = y_values[i];
629+
yp[i] = yp_values[i];
605630
}
631+
model.y().setDataUpdated();
632+
model.yp().setDataUpdated();
606633

607634
model.updateTime(0.0, 1.0);
608635
model.evaluateResidual();

0 commit comments

Comments
 (0)