Inequality Jacobian indexes are correct (for lines)

pnnl · Oct 31, 2023 · 136c505 · 136c505
1 parent 365eae7
commit 136c505
Show file tree

Hide file tree

Showing 4 changed files with 188 additions and 83 deletions.
diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp
@@ -235,6 +235,7 @@ int LINEParamsRajaHiop::copy(OPFLOW opflow) {
     resmgr.copy(gineqidx_dev_, gineqidx);
     resmgr.copy(gbineqidx_dev_, gbineqidx);
     resmgr.copy(linelimidx_dev_, linelimidx);
+    resmgr.copy(jac_ieq_idx_dev_, jac_ieq_idx);
   }
 #else
   Gff_dev_ = Gff;
@@ -258,6 +259,7 @@ int LINEParamsRajaHiop::copy(OPFLOW opflow) {
     gineqidx_dev_ = gineqidx;
     gbineqidx_dev_ = gbineqidx;
     linelimidx_dev_ = linelimidx;
+    jac_ieq_idx_dev_ = jac_ieq_idx;
   }
 #endif
   return 0;
@@ -289,6 +291,7 @@ int LINEParamsRajaHiop::destroy(OPFLOW opflow) {
     h_allocator_.deallocate(gineqidx);
     h_allocator_.deallocate(gbineqidx);
     h_allocator_.deallocate(linelimidx);
+    h_allocator_.deallocate(jac_ieq_idx);
   }
 
 #ifdef EXAGO_ENABLE_GPU
@@ -317,6 +320,7 @@ int LINEParamsRajaHiop::destroy(OPFLOW opflow) {
     d_allocator_.deallocate(gineqidx_dev_);
     d_allocator_.deallocate(gbineqidx_dev_);
     d_allocator_.deallocate(linelimidx_dev_);
+    d_allocator_.deallocate(jac_ieq_idx_dev_);
   }
 #endif
 
@@ -369,6 +373,7 @@ int LINEParamsRajaHiop::allocate(OPFLOW opflow) {
     linelimidx = paramAlloc<int>(h_allocator_, nlinelim);
     gineqidx = paramAlloc<int>(h_allocator_, nlinelim);
     gbineqidx = paramAlloc<int>(h_allocator_, nlinelim);
+    jac_ieq_idx = paramAlloc<int>(h_allocator_, nlinelim);
   }
 
   PetscInt j = 0;
@@ -416,6 +421,7 @@ int LINEParamsRajaHiop::allocate(OPFLOW opflow) {
       gbineqidx[j] = opflow->nconeq + line->startineqloc;
       gineqidx[j] = line->startineqloc;
       linelimidx[j] = linei;
+      jac_ieq_idx[j] = 0;
       j++;
     }
 
@@ -450,6 +456,7 @@ int LINEParamsRajaHiop::allocate(OPFLOW opflow) {
     gineqidx_dev_ = paramAlloc<int>(d_allocator_, nlinelim);
     gbineqidx_dev_ = paramAlloc<int>(d_allocator_, nlinelim);
     linelimidx_dev_ = paramAlloc<int>(d_allocator_, nlinelim);
+    jac_ieq_idx_dev_ = paramAlloc<int>(d_allocator_, nlinelim);
   }
 #endif
   return 0;

diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.h b/src/opflow/model/power_bal_hiop/paramsrajahiop.h
@@ -200,6 +200,7 @@ struct LINEParamsRajaHiop {
   int *bust_idx;   /* To bus index */
   int *jacf_idx;   /* Location number in the sparse Jacobian (from) */
   int *jact_idx;   /* Location number in the sparse Jacobian (to) */
+  int *jac_ieq_idx;/* Location number in sparse inequality Jacobian */
 
   // Device data
   double *Gff_dev_;    /* From side self conductance */
@@ -228,6 +229,7 @@ struct LINEParamsRajaHiop {
   int *bust_idx_dev_;   /* To bus index */
   int *jacf_idx_dev_;   /* Location number in the sparse Jacobian (from) */
   int *jact_idx_dev_;   /* Location number in the sparse Jacobian (to) */
+  int *jac_ieq_idx_dev_;/* Location number in sparse inequality Jacobian */
 
   int allocate(OPFLOW);
   int destroy(OPFLOW);

diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp
@@ -310,19 +310,25 @@ PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) {
     nnz_eqjac += 4;
   }
 
+  // if there are lines, non-zeros were over counted
+  if (ps->nline > 0) {
+    nnz_eqjac -= 8;
+  }
+
   std::cout << "Equality Jacobian nonzero count: " << nnz_eqjac << std::endl;
 
   if (opflow->has_gensetpoint) {
-    for (int ibus = 0; ibus < ps->nbus; ++ibus) {
+    for (int ibus = 0, igen = 0; ibus < ps->nbus; ++ibus) {
       PSBUS bus = &(ps->bus[ibus]);
-      for (int igen = 0; igen < bus->ngen; ++igen) {
+      for (int bgen = 0; bgen < bus->ngen; ++bgen) {
         PSGEN gen;
-        ierr = PSBUSGetGen(bus, igen, &gen);
+        ierr = PSBUSGetGen(bus, bgen, &gen);
         CHKERRQ(ierr);
         if (!gen->status)
           continue;
-        genparams->ineqjacspgen_idx[igen] = nnz_ineqjac;
+        genparams->ineqjacspgen_idx[igen] = nnz_eqjac + nnz_ineqjac;
         nnz_ineqjac += 6;
+        igen++;
       }
     }
   }
@@ -331,9 +337,9 @@ PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) {
     for (int ibus = 0; ibus < ps->nbus; ++ibus) {
       PSBUS bus = &(ps->bus[ibus]);
       if (bus->ide == PV_BUS || bus->ide == REF_BUS) {
-        for (int igen = 0; igen < bus->ngen; ++igen) {
+        for (int bgen = 0; bgen < bus->ngen; ++bgen) {
         PSGEN gen;
-        ierr = PSBUSGetGen(bus, igen, &gen);
+        ierr = PSBUSGetGen(bus, bgen, &gen);
         CHKERRQ(ierr);
         if (!gen->status)
           continue;
@@ -346,6 +352,7 @@ PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) {
   if (!opflow->ignore_lineflow_constraints) {
     for (int iline = 0; iline < opflow->nlinesmon; ++iline) {
       // PSLINE line = &ps->line[opflow->linesmon[iline]];
+      lineparams->jac_ieq_idx[iline] = nnz_eqjac + nnz_ineqjac;
       nnz_ineqjac += 8;
     }
   }