diff --git a/include/common.h b/include/common.h index 511f79f5..357a2e58 100644 --- a/include/common.h +++ b/include/common.h @@ -7,6 +7,7 @@ #define COMMON_H #include +#include typedef enum { MATPOWER, CSV, JSON, MINIMAL } OutputFormat; diff --git a/include/opflow.h b/include/opflow.h index 7c1c580c..e1b1ae45 100644 --- a/include/opflow.h +++ b/include/opflow.h @@ -16,11 +16,13 @@ #define OPFLOWMODEL_PBPOLHIOP "POWER_BALANCE_HIOP" #define OPFLOWMODEL_PBPOLRAJAHIOP "PBPOLRAJAHIOP" #define OPFLOWMODEL_DCOPF "DCOPF" +#define OPFLOWMODEL_PBPOLRAJAHIOPSPARSE "PBPOLRAJAHIOPSPARSE" /* Solvers */ #define OPFLOWSOLVER_IPOPT "IPOPT" #define OPFLOWSOLVER_HIOP "HIOP" #define OPFLOWSOLVER_HIOPSPARSE "HIOPSPARSE" +#define OPFLOWSOLVER_HIOPSPARSEGPU "HIOPSPARSEGPU" typedef struct _p_OPFLOW *OPFLOW; diff --git a/src/opflow/CMakeLists.txt b/src/opflow/CMakeLists.txt index b6fe1592..b4dd9e0b 100644 --- a/src/opflow/CMakeLists.txt +++ b/src/opflow/CMakeLists.txt @@ -11,13 +11,22 @@ set(OPFLOW_FORM_SRC # Build RAJA/Umpire modules if RAJA and Umpire are enabled if(EXAGO_ENABLE_RAJA) - set(OPFLOW_FORM_SRC ${OPFLOW_FORM_SRC} model/power_bal_hiop/pbpolrajahiop.cpp - model/power_bal_hiop/pbpolrajahiopkernels.cpp + set(OPFLOW_FORM_SRC + ${OPFLOW_FORM_SRC} model/power_bal_hiop/pbpolrajahiop.cpp + model/power_bal_hiop/pbpolrajahiopkernels.cpp + model/power_bal_hiop/paramsrajahiop.cpp ) + if(EXAGO_ENABLE_HIOP_SPARSE) + set(OPFLOW_FORM_SRC + ${OPFLOW_FORM_SRC} model/power_bal_hiop/pbpolrajahiopsparse.cpp + model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp + ) + endif() endif() -set(OPFLOW_SOLVER_SRC solver/ipopt/opflow_ipopt.cpp solver/hiop/opflow_hiop.cpp - solver/hiop/opflow_hiopsparse.cpp +set(OPFLOW_SOLVER_SRC + solver/ipopt/opflow_ipopt.cpp solver/hiop/opflow_hiop.cpp + solver/hiop/opflow_hiopsparse.cpp solver/hiop/opflow_hiopsparsegpu.cpp ) set(OPFLOW_SRC @@ -29,7 +38,9 @@ set_source_files_properties(${OPFLOW_SRC} PROPERTIES LANGUAGE CXX) if(EXAGO_ENABLE_RAJA AND EXAGO_ENABLE_CUDA) set_source_files_properties( - model/power_bal_hiop/pbpolrajahiopkernels.cpp PROPERTIES LANGUAGE CUDA + model/power_bal_hiop/pbpolrajahiopkernels.cpp + model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp PROPERTIES LANGUAGE + CUDA ) endif() diff --git a/src/opflow/interface/opflow.cpp b/src/opflow/interface/opflow.cpp index c3f97b78..2c6787c2 100644 --- a/src/opflow/interface/opflow.cpp +++ b/src/opflow/interface/opflow.cpp @@ -3008,6 +3008,11 @@ PetscErrorCode OPFLOWCheckModelSolverCompatibility(OPFLOW opflow) { #else PetscBool rajahiop_pbpol = PETSC_FALSE; #endif // RAJA + if (hiop && !(hiop_pbpol || rajahiop_pbpol)) { + SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, + "OPFLOW solver HIOP incompatible with model %s", + (opflow->modelname).c_str()); + } #if defined(EXAGO_ENABLE_HIOP_SPARSE) PetscBool hiop_sparse; PetscBool hiop_sparse_pbpol; @@ -3022,14 +3027,22 @@ PetscErrorCode OPFLOWCheckModelSolverCompatibility(OPFLOW opflow) { "OPFLOW solver HIOPSPARSE incompatible with model %s", (opflow->modelname).c_str()); } -#else - PetscBool hiop_sparse = PETSC_FALSE; -#endif // HIOP_SPARSE - if ((hiop && !hiop_sparse) && !(hiop_pbpol || rajahiop_pbpol)) { +#if defined(EXAGO_ENABLE_RAJA) + PetscBool hiop_sparsegpu; + PetscBool pbpolrajahiopsparse; + + hiop_sparsegpu = + static_cast(opflow->solvername == OPFLOWSOLVER_HIOPSPARSEGPU); + pbpolrajahiopsparse = static_cast(opflow->modelname == + OPFLOWMODEL_PBPOLRAJAHIOPSPARSE); + + if (hiop_sparsegpu && !pbpolrajahiopsparse) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, - "OPFLOW solver HIOP incompatible with model %s", + "OPFLOW solver HIOPSPARSE incompatible with model %s", (opflow->modelname).c_str()); } +#endif // EXAGO_ENABLE_RAJA +#endif // HIOP_SPARSE /* FIXME: The PBPOLRAJAHIOP has trouble with individual load shedding, so avoid using when load shedding is enabled */ if (opflow->include_loadloss_variables) { diff --git a/src/opflow/interface/opflowregi.cpp b/src/opflow/interface/opflowregi.cpp index a9b16a9d..f72698b5 100644 --- a/src/opflow/interface/opflowregi.cpp +++ b/src/opflow/interface/opflowregi.cpp @@ -76,6 +76,7 @@ extern PetscErrorCode OPFLOWModelCreate_PBPOLHIOP(OPFLOW); #if defined(EXAGO_ENABLE_RAJA) extern PetscErrorCode OPFLOWModelCreate_PBPOLRAJAHIOP(OPFLOW); +extern PetscErrorCode OPFLOWModelCreate_PBPOLRAJAHIOPSPARSE(OPFLOW); #endif /* @@ -114,6 +115,9 @@ PetscErrorCode OPFLOWModelRegisterAll(OPFLOW opflow) { ierr = OPFLOWModelRegister(opflow, OPFLOWMODEL_PBPOLRAJAHIOP, OPFLOWModelCreate_PBPOLRAJAHIOP); CHKERRQ(ierr); + ierr = OPFLOWModelRegister(opflow, OPFLOWMODEL_PBPOLRAJAHIOPSPARSE, + OPFLOWModelCreate_PBPOLRAJAHIOPSPARSE); + CHKERRQ(ierr); #endif opflow->OPFLOWModelRegisterAllCalled = PETSC_TRUE; @@ -128,6 +132,9 @@ extern PetscErrorCode OPFLOWSolverCreate_IPOPT(OPFLOW); extern PetscErrorCode OPFLOWSolverCreate_HIOP(OPFLOW); #if defined(EXAGO_ENABLE_HIOP_SPARSE) extern PetscErrorCode OPFLOWSolverCreate_HIOPSPARSE(OPFLOW); +#if defined(EXAGO_ENABLE_RAJA) +extern PetscErrorCode OPFLOWSolverCreate_HIOPSPARSEGPU(OPFLOW); +#endif #endif #endif @@ -153,6 +160,11 @@ PetscErrorCode OPFLOWSolverRegisterAll(OPFLOW opflow) { ierr = OPFLOWSolverRegister(opflow, OPFLOWSOLVER_HIOPSPARSE, OPFLOWSolverCreate_HIOPSPARSE); CHKERRQ(ierr); +#if defined(EXAGO_ENABLE_RAJA) + ierr = OPFLOWSolverRegister(opflow, OPFLOWSOLVER_HIOPSPARSEGPU, + OPFLOWSolverCreate_HIOPSPARSEGPU); + CHKERRQ(ierr); +#endif #endif #endif opflow->OPFLOWSolverRegisterAllCalled = PETSC_TRUE; diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp new file mode 100644 index 00000000..27e2cfa0 --- /dev/null +++ b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp @@ -0,0 +1,800 @@ + +#include "paramsrajahiop.h" +#include + +template T *paramAlloc(umpire::Allocator &a, I N) { + return static_cast(a.allocate(N * sizeof(T))); +} + +/* Functions to create and destroy data arrays for different + component classes +*/ +int BUSParamsRajaHiop::destroy(OPFLOW opflow) { + h_allocator_.deallocate(isref); + h_allocator_.deallocate(isisolated); + h_allocator_.deallocate(ispvpq); + h_allocator_.deallocate(vmin); + h_allocator_.deallocate(vmax); + h_allocator_.deallocate(va); + h_allocator_.deallocate(vm); + h_allocator_.deallocate(gl); + h_allocator_.deallocate(bl); + h_allocator_.deallocate(xidx); + h_allocator_.deallocate(gidx); + if (opflow->include_powerimbalance_variables) { + h_allocator_.deallocate(xidxpimb); + h_allocator_.deallocate(powerimbalance_penalty); + h_allocator_.deallocate(jacsp_idx); + h_allocator_.deallocate(jacsq_idx); + } + +#ifdef EXAGO_ENABLE_GPU + d_allocator_.deallocate(isref_dev_); + d_allocator_.deallocate(isisolated_dev_); + d_allocator_.deallocate(ispvpq_dev_); + d_allocator_.deallocate(vmin_dev_); + d_allocator_.deallocate(vmax_dev_); + d_allocator_.deallocate(va_dev_); + d_allocator_.deallocate(vm_dev_); + d_allocator_.deallocate(gl_dev_); + d_allocator_.deallocate(bl_dev_); + d_allocator_.deallocate(xidx_dev_); + d_allocator_.deallocate(gidx_dev_); + if (opflow->include_powerimbalance_variables) { + d_allocator_.deallocate(xidxpimb_dev_); + d_allocator_.deallocate(powerimbalance_penalty_dev_); + d_allocator_.deallocate(jacsp_idx_dev_); + d_allocator_.deallocate(jacsq_idx_dev_); + } +#endif + + return 0; +} + +/* Copy data from host to device */ +int BUSParamsRajaHiop::copy(OPFLOW opflow) { + /* Allocate the arrays */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + +#ifdef EXAGO_ENABLE_GPU + d_allocator_ = resmgr.getAllocator("DEVICE"); + + // Copy host data from the host to the device + resmgr.copy(isref_dev_, isref); + resmgr.copy(isisolated_dev_, isisolated); + resmgr.copy(ispvpq_dev_, ispvpq); + + resmgr.copy(vmin_dev_, vmin); + resmgr.copy(vmax_dev_, vmax); + resmgr.copy(va_dev_, va); + resmgr.copy(vm_dev_, vm); + resmgr.copy(gl_dev_, gl); + resmgr.copy(bl_dev_, bl); + + resmgr.copy(xidx_dev_, xidx); + resmgr.copy(gidx_dev_, gidx); + if (opflow->include_powerimbalance_variables) { + resmgr.copy(xidxpimb_dev_, xidxpimb); + resmgr.copy(jacsp_idx_dev_, jacsp_idx); + resmgr.copy(jacsq_idx_dev_, jacsq_idx); + resmgr.copy(powerimbalance_penalty_dev_, powerimbalance_penalty); + } +#else + isref_dev_ = isref; + isisolated_dev_ = isisolated; + ispvpq_dev_ = ispvpq; + vmin_dev_ = vmin; + vmax_dev_ = vmax; + va_dev_ = va; + vm_dev_ = vm; + gl_dev_ = gl; + bl_dev_ = bl; + xidx_dev_ = xidx; + xidxpimb_dev_ = xidxpimb; + gidx_dev_ = gidx; + jacsp_idx_dev_ = jacsp_idx; + jacsq_idx_dev_ = jacsq_idx; + powerimbalance_penalty_dev_ = powerimbalance_penalty; +#endif + return 0; +} + +/* Create data for buses that is used in different computations */ +int BUSParamsRajaHiop::allocate(OPFLOW opflow) { + PS ps = opflow->ps; + PetscInt loc; + PSBUS bus; + + nbus = ps->nbus; + /* Allocate the arrays */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + + // Allocate data on the host + isref = paramAlloc(h_allocator_, nbus); + isisolated = paramAlloc(h_allocator_, nbus); + ispvpq = paramAlloc(h_allocator_, nbus); + + vmin = paramAlloc(h_allocator_, nbus); + vmax = paramAlloc(h_allocator_, nbus); + va = paramAlloc(h_allocator_, nbus); + vm = paramAlloc(h_allocator_, nbus); + gl = paramAlloc(h_allocator_, nbus); + bl = paramAlloc(h_allocator_, nbus); + + xidx = paramAlloc(h_allocator_, nbus); + gidx = paramAlloc(h_allocator_, nbus); + + if (opflow->include_powerimbalance_variables) { + xidxpimb = paramAlloc(h_allocator_, nbus); + powerimbalance_penalty = paramAlloc(h_allocator_, nbus); + jacsp_idx = paramAlloc(h_allocator_, nbus); + jacsq_idx = paramAlloc(h_allocator_, nbus); + } + + /* Memzero arrays */ + resmgr.memset(isref, 0, nbus * sizeof(int)); + resmgr.memset(ispvpq, 0, nbus * sizeof(int)); + resmgr.memset(isisolated, 0, nbus * sizeof(int)); + + for (int i = 0; i < nbus; i++) { + bus = &ps->bus[i]; + loc = bus->startxVloc; + + xidx[i] = opflow->idxn2sd_map[loc]; + gidx[i] = bus->starteqloc; + + if (bus->ide == REF_BUS) + isref[i] = 1; + else if (bus->ide == ISOLATED_BUS) + isisolated[i] = 1; + else + ispvpq[i] = 1; + + if (opflow->genbusvoltagetype == FIXED_AT_SETPOINT) { + if (bus->ide == REF_BUS || bus->ide == PV_BUS) { + /* Hold voltage at reference and PV buses */ + vmin[i] = bus->vm; + vmax[i] = bus->vm; + } else { + vmin[i] = bus->Vmin; + vmax[i] = bus->Vmax; + } + } else { + vmin[i] = bus->Vmin; + vmax[i] = bus->Vmax; + } + vm[i] = bus->vm; + va[i] = bus->va; + gl[i] = bus->gl; + bl[i] = bus->bl; + + if (opflow->include_powerimbalance_variables) { + loc = bus->startxpimbloc; + xidxpimb[i] = opflow->idxn2sd_map[loc]; + powerimbalance_penalty[i] = opflow->powerimbalance_penalty; + } + } + +#ifdef EXAGO_ENABLE_GPU + d_allocator_ = resmgr.getAllocator("DEVICE"); + // Allocate data on the device + isref_dev_ = paramAlloc(d_allocator_, nbus); + isisolated_dev_ = paramAlloc(d_allocator_, nbus); + ispvpq_dev_ = paramAlloc(d_allocator_, nbus); + + vmin_dev_ = paramAlloc(d_allocator_, nbus); + vmax_dev_ = paramAlloc(d_allocator_, nbus); + va_dev_ = paramAlloc(d_allocator_, nbus); + vm_dev_ = paramAlloc(d_allocator_, nbus); + gl_dev_ = paramAlloc(d_allocator_, nbus); + bl_dev_ = paramAlloc(d_allocator_, nbus); + + xidx_dev_ = paramAlloc(d_allocator_, nbus); + gidx_dev_ = paramAlloc(d_allocator_, nbus); + + if (opflow->include_powerimbalance_variables) { + xidxpimb_dev_ = paramAlloc(d_allocator_, nbus); + powerimbalance_penalty_dev_ = paramAlloc(d_allocator_, nbus); + jacsp_idx_dev_ = paramAlloc(d_allocator_, nbus); + jacsq_idx_dev_ = paramAlloc(d_allocator_, nbus); + } +#endif + return 0; +} + +int LINEParamsRajaHiop::copy(OPFLOW opflow) { + /* Allocate the arrays */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + +#ifdef EXAGO_ENABLE_GPU + // Copy data from the host to the device + resmgr.copy(Gff_dev_, Gff); + resmgr.copy(Bff_dev_, Bff); + resmgr.copy(Gft_dev_, Gft); + resmgr.copy(Bft_dev_, Bft); + resmgr.copy(Gtf_dev_, Gtf); + resmgr.copy(Btf_dev_, Btf); + resmgr.copy(Gtt_dev_, Gtt); + resmgr.copy(Btt_dev_, Btt); + resmgr.copy(rateA_dev_, rateA); + + resmgr.copy(xidxf_dev_, xidxf); + resmgr.copy(xidxt_dev_, xidxt); + + resmgr.copy(geqidxf_dev_, geqidxf); + resmgr.copy(geqidxt_dev_, geqidxt); + + if (opflow->nlinesmon) { + resmgr.copy(gineqidx_dev_, gineqidx); + resmgr.copy(gbineqidx_dev_, gbineqidx); + resmgr.copy(linelimidx_dev_, linelimidx); + } +#else + Gff_dev_ = Gff; + Bff_dev_ = Bff; + Gft_dev_ = Gft; + Bft_dev_ = Bft; + Gtf_dev_ = Gtf; + Btf_dev_ = Btf; + Gtt_dev_ = Gtt; + Btt_dev_ = Btt; + rateA_dev_ = rateA; + xidxf_dev_ = xidxf; + xidxt_dev_ = xidxt; + geqidxf_dev_ = geqidxf; + geqidxt_dev_ = geqidxt; + if (opflow->nlinesmon) { + gineqidx_dev_ = gineqidx; + gbineqidx_dev_ = gbineqidx; + linelimidx_dev_ = linelimidx; + } +#endif + return 0; +} + +int LINEParamsRajaHiop::destroy(OPFLOW opflow) { + // Destroy parameter arrays on the host + h_allocator_.deallocate(Gff); + h_allocator_.deallocate(Bff); + h_allocator_.deallocate(Gft); + h_allocator_.deallocate(Bft); + h_allocator_.deallocate(Gtf); + h_allocator_.deallocate(Btf); + h_allocator_.deallocate(Gtt); + h_allocator_.deallocate(Btt); + h_allocator_.deallocate(rateA); + + h_allocator_.deallocate(xidxf); + h_allocator_.deallocate(xidxt); + + h_allocator_.deallocate(geqidxf); + h_allocator_.deallocate(geqidxt); + + if (opflow->nlinesmon) { + h_allocator_.deallocate(gineqidx); + h_allocator_.deallocate(gbineqidx); + h_allocator_.deallocate(linelimidx); + } + +#ifdef EXAGO_ENABLE_GPU + // Destroy parameter arrays on the device + d_allocator_.deallocate(Gff_dev_); + d_allocator_.deallocate(Bff_dev_); + d_allocator_.deallocate(Gft_dev_); + d_allocator_.deallocate(Bft_dev_); + d_allocator_.deallocate(Gtf_dev_); + d_allocator_.deallocate(Btf_dev_); + d_allocator_.deallocate(Gtt_dev_); + d_allocator_.deallocate(Btt_dev_); + d_allocator_.deallocate(rateA_dev_); + + d_allocator_.deallocate(xidxf_dev_); + d_allocator_.deallocate(xidxt_dev_); + + d_allocator_.deallocate(geqidxf_dev_); + d_allocator_.deallocate(geqidxt_dev_); + + if (opflow->nlinesmon) { + d_allocator_.deallocate(gineqidx_dev_); + d_allocator_.deallocate(gbineqidx_dev_); + d_allocator_.deallocate(linelimidx_dev_); + } +#endif + + return 0; +} + +/* Create data for lines that is used in different computations */ +int LINEParamsRajaHiop::allocate(OPFLOW opflow) { + PS ps = opflow->ps; + PetscInt linei = 0; + PSLINE line; + PetscInt i; + PetscErrorCode ierr; + const PSBUS *connbuses; + PSBUS busf, bust; + + PetscFunctionBegin; + ierr = PSGetNumActiveLines(ps, &nlineON, NULL); + CHKERRQ(ierr); + + nlinelim = opflow->nlinesmon; + + /* Allocate data arrays */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + + // Allocate data on the host + Gff = paramAlloc(h_allocator_, nlineON); + Bff = paramAlloc(h_allocator_, nlineON); + Gft = paramAlloc(h_allocator_, nlineON); + Bft = paramAlloc(h_allocator_, nlineON); + Gtf = paramAlloc(h_allocator_, nlineON); + Btf = paramAlloc(h_allocator_, nlineON); + Gtt = paramAlloc(h_allocator_, nlineON); + Btt = paramAlloc(h_allocator_, nlineON); + rateA = paramAlloc(h_allocator_, nlineON); + + xidxf = paramAlloc(h_allocator_, nlineON); + xidxt = paramAlloc(h_allocator_, nlineON); + + geqidxf = paramAlloc(h_allocator_, nlineON); + geqidxt = paramAlloc(h_allocator_, nlineON); + + if (opflow->nlinesmon) { + linelimidx = paramAlloc(h_allocator_, nlinelim); + gineqidx = paramAlloc(h_allocator_, nlinelim); + gbineqidx = paramAlloc(h_allocator_, nlinelim); + } + + PetscInt j = 0; + /* Populate arrays */ + for (i = 0; i < ps->nline; i++) { + line = &ps->line[i]; + + if (!line->status) + continue; + + Gff[linei] = line->yff[0]; + Bff[linei] = line->yff[1]; + Gft[linei] = line->yft[0]; + Bft[linei] = line->yft[1]; + Gtf[linei] = line->ytf[0]; + Btf[linei] = line->ytf[1]; + Gtt[linei] = line->ytt[0]; + Btt[linei] = line->ytt[1]; + rateA[linei] = line->rateA; + + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; + + int xidxfi, xidxti; + xidxfi = busf->startxVloc; + xidxti = bust->startxVloc; + + xidxf[linei] = opflow->idxn2sd_map[xidxfi]; + xidxt[linei] = opflow->idxn2sd_map[xidxti]; + + /* + Each bus has two equality (balance) constraints, hence the use of + coefficient 2 to map the location of the equality constraint for the bus + */ + geqidxf[linei] = busf->starteqloc; + geqidxt[linei] = bust->starteqloc; + + if (j < opflow->nlinesmon && opflow->linesmon[j] == i) { + gbineqidx[j] = opflow->nconeq + line->startineqloc; + gineqidx[j] = line->startineqloc; + linelimidx[j] = linei; + j++; + } + + linei++; + } + +#ifdef EXAGO_ENABLE_GPU + d_allocator_ = resmgr.getAllocator("DEVICE"); + // Allocate data on the device + Gff_dev_ = paramAlloc(d_allocator_, nlineON); + Bff_dev_ = paramAlloc(d_allocator_, nlineON); + Gft_dev_ = paramAlloc(d_allocator_, nlineON); + Bft_dev_ = paramAlloc(d_allocator_, nlineON); + Gtf_dev_ = paramAlloc(d_allocator_, nlineON); + Btf_dev_ = paramAlloc(d_allocator_, nlineON); + Gtt_dev_ = paramAlloc(d_allocator_, nlineON); + Btt_dev_ = paramAlloc(d_allocator_, nlineON); + rateA_dev_ = paramAlloc(d_allocator_, nlineON); + + xidxf_dev_ = paramAlloc(d_allocator_, nlineON); + xidxt_dev_ = paramAlloc(d_allocator_, nlineON); + + geqidxf_dev_ = paramAlloc(d_allocator_, nlineON); + geqidxt_dev_ = paramAlloc(d_allocator_, nlineON); + + if (opflow->nconineq) { + gineqidx_dev_ = paramAlloc(d_allocator_, nlinelim); + gbineqidx_dev_ = paramAlloc(d_allocator_, nlinelim); + linelimidx_dev_ = paramAlloc(d_allocator_, nlinelim); + } +#endif + return 0; +} + +int LOADParamsRajaHiop::copy(OPFLOW opflow) { + /* Allocate the arrays */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + +#ifdef EXAGO_ENABLE_GPU + // Copy data from host to device + resmgr.copy(pl_dev_, pl); + resmgr.copy(ql_dev_, ql); + resmgr.copy(xidx_dev_, xidx); + resmgr.copy(gidx_dev_, gidx); + if (opflow->include_loadloss_variables) { + resmgr.copy(jacsp_idx_dev_, jacsp_idx); + resmgr.copy(jacsq_idx_dev_, jacsq_idx); + resmgr.copy(hesssp_idx_dev_, hesssp_idx); + resmgr.copy(loadloss_penalty_dev_, loadloss_penalty); + } +#else + pl_dev_ = pl; + ql_dev_ = ql; + xidx_dev_ = xidx; + gidx_dev_ = gidx; + jacsp_idx_dev_ = jacsp_idx; + jacsq_idx_dev_ = jacsq_idx; + hesssp_idx_dev_ = hesssp_idx; + loadloss_penalty_dev_ = loadloss_penalty; +#endif + + return 0; +} + +int LOADParamsRajaHiop::destroy(OPFLOW opflow) { + h_allocator_.deallocate(pl); + h_allocator_.deallocate(ql); + h_allocator_.deallocate(xidx); + h_allocator_.deallocate(gidx); + if (opflow->include_loadloss_variables) { + h_allocator_.deallocate(loadloss_penalty); + h_allocator_.deallocate(jacsp_idx); + h_allocator_.deallocate(jacsq_idx); + h_allocator_.deallocate(hesssp_idx); + } +#ifdef EXAGO_ENABLE_GPU + d_allocator_.deallocate(pl_dev_); + d_allocator_.deallocate(ql_dev_); + d_allocator_.deallocate(xidx_dev_); + d_allocator_.deallocate(gidx_dev_); + if (opflow->include_loadloss_variables) { + d_allocator_.deallocate(loadloss_penalty_dev_); + d_allocator_.deallocate(jacsp_idx_dev_); + d_allocator_.deallocate(jacsq_idx_dev_); + d_allocator_.deallocate(hesssp_idx_dev_); + } +#endif + return 0; +} + +/* Create data for loads that is used in different computations */ +int LOADParamsRajaHiop::allocate(OPFLOW opflow) { + PS ps = opflow->ps; + PetscInt loc, loadi = 0; + PSLOAD load; + PSBUS bus; + PetscInt i, j; + PetscErrorCode ierr; + + /* Get the number of active generators (STATUS ON) */ + ierr = PSGetNumLoads(ps, &nload, NULL); + CHKERRQ(ierr); + + /* Allocate arrays */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + + // Allocate data on the host + pl = paramAlloc(h_allocator_, nload); + ql = paramAlloc(h_allocator_, nload); + xidx = paramAlloc(h_allocator_, nload); + gidx = paramAlloc(h_allocator_, nload); + if (opflow->include_loadloss_variables) { + loadloss_penalty = paramAlloc(h_allocator_, nload); + jacsp_idx = paramAlloc(h_allocator_, nload); + jacsq_idx = paramAlloc(h_allocator_, nload); + hesssp_idx = paramAlloc(h_allocator_, nload); + } + /* Insert data in loadparams */ + for (i = 0; i < ps->nbus; i++) { + bus = &ps->bus[i]; + ierr = PSBUSGetVariableLocation(bus, &loc); + CHKERRQ(ierr); + for (j = 0; j < bus->nload; j++) { + ierr = PSBUSGetLoad(bus, j, &load); + CHKERRQ(ierr); + pl[loadi] = load->pl; + ql[loadi] = load->ql; + if (opflow->include_loadloss_variables) { + loc = load->startxloadlossloc; + loadloss_penalty[loadi] = load->loss_cost; + xidx[loadi] = opflow->idxn2sd_map[loc]; + } + gidx[loadi] = bus->starteqloc; + loadi++; + } + } +#ifdef EXAGO_ENABLE_GPU + d_allocator_ = resmgr.getAllocator("DEVICE"); + // Allocate data on the device + pl_dev_ = paramAlloc(d_allocator_, nload); + ql_dev_ = paramAlloc(d_allocator_, nload); + xidx_dev_ = paramAlloc(d_allocator_, nload); + gidx_dev_ = paramAlloc(d_allocator_, nload); + if (opflow->include_loadloss_variables) { + loadloss_penalty_dev_ = paramAlloc(d_allocator_, nload); + jacsp_idx_dev_ = paramAlloc(d_allocator_, nload); + jacsq_idx_dev_ = paramAlloc(d_allocator_, nload); + hesssp_idx_dev_ = paramAlloc(d_allocator_, nload); + } + +#endif + + return (0); +} + +int GENParamsRajaHiop::destroy(OPFLOW opflow) { + // Free arrays on the host + h_allocator_.deallocate(cost_alpha); + h_allocator_.deallocate(cost_beta); + h_allocator_.deallocate(cost_gamma); + h_allocator_.deallocate(pt); + h_allocator_.deallocate(pb); + h_allocator_.deallocate(qt); + h_allocator_.deallocate(qb); + h_allocator_.deallocate(isrenewable); + h_allocator_.deallocate(xidx); + h_allocator_.deallocate(gidxbus); + h_allocator_.deallocate(eqjacspbus_idx); + h_allocator_.deallocate(eqjacsqbus_idx); + h_allocator_.deallocate(hesssp_idx); + if (opflow->has_gensetpoint) { + h_allocator_.deallocate(geqidxgen); + h_allocator_.deallocate(gineqidxgen); + h_allocator_.deallocate(gbineqidxgen); + h_allocator_.deallocate(pgs); + h_allocator_.deallocate(eqjacspgen_idx); + h_allocator_.deallocate(ineqjacspgen_idx); + } +#ifdef EXAGO_ENABLE_GPU + // Free arrays on the device + d_allocator_.deallocate(cost_alpha_dev_); + d_allocator_.deallocate(cost_beta_dev_); + d_allocator_.deallocate(cost_gamma_dev_); + d_allocator_.deallocate(pt_dev_); + d_allocator_.deallocate(pb_dev_); + d_allocator_.deallocate(qt_dev_); + d_allocator_.deallocate(qb_dev_); + d_allocator_.deallocate(isrenewable_dev_); + d_allocator_.deallocate(xidx_dev_); + d_allocator_.deallocate(gidxbus_dev_); + d_allocator_.deallocate(eqjacspbus_idx_dev_); + d_allocator_.deallocate(eqjacsqbus_idx_dev_); + d_allocator_.deallocate(hesssp_idx_dev_); + if (opflow->has_gensetpoint) { + d_allocator_.deallocate(geqidxgen_dev_); + d_allocator_.deallocate(gineqidxgen_dev_); + d_allocator_.deallocate(gbineqidxgen_dev_); + d_allocator_.deallocate(pgs_dev_); + d_allocator_.deallocate(eqjacspgen_idx_dev_); + d_allocator_.deallocate(ineqjacspgen_idx_dev_); + } + +#endif + return 0; +} + +int GENParamsRajaHiop::copy(OPFLOW opflow) { + /* Allocate arrays on the host */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + +#ifdef EXAGO_ENABLE_GPU + // Copy host data to the device + resmgr.copy(cost_alpha_dev_, cost_alpha); + resmgr.copy(cost_beta_dev_, cost_beta); + resmgr.copy(cost_gamma_dev_, cost_gamma); + + resmgr.copy(pt_dev_, pt); + resmgr.copy(pb_dev_, pb); + resmgr.copy(qt_dev_, qt); + resmgr.copy(qb_dev_, qb); + resmgr.copy(isrenewable_dev_, isrenewable); + + resmgr.copy(xidx_dev_, xidx); + resmgr.copy(gidxbus_dev_, gidxbus); + + resmgr.copy(eqjacspbus_idx_dev_, eqjacspbus_idx); + resmgr.copy(eqjacsqbus_idx_dev_, eqjacsqbus_idx); + resmgr.copy(hesssp_idx_dev_, hesssp_idx); + if (opflow->has_gensetpoint) { + resmgr.copy(geqidxgen_dev_, geqidxgen); + resmgr.copy(gineqidxgen_dev_, gineqidxgen); + resmgr.copy(gbineqidxgen_dev_, gbineqidxgen); + resmgr.copy(eqjacspgen_idx_dev_, eqjacspgen_idx); + resmgr.copy(ineqjacspgen_idx_dev_, ineqjacspgen_idx); + resmgr.copy(pgs_dev_, pgs); + } +#else + cost_alpha_dev_ = cost_alpha; + cost_beta_dev_ = cost_beta; + cost_gamma_dev_ = cost_gamma; + pt_dev_ = pt; + pb_dev_ = pb; + qt_dev_ = qt; + qb_dev_ = qb; + isrenewable_dev_ = isrenewable; + xidx_dev_ = xidx; + gidxbus_dev_ = gidxbus; + eqjacspbus_idx_dev_ = eqjacspbus_idx; + eqjacsqbus_idx_dev_ = eqjacsqbus_idx; + hesssp_idx_dev_ = hesssp_idx; + geqidxgen_dev_ = geqidxgen; + gineqidxgen_dev_ = gineqidxgen; + gbineqidxgen_dev_ = gbineqidxgen; + eqjacspgen_idx_dev_ = eqjacspgen_idx; + ineqjacspgen_idx_dev_ = ineqjacspgen_idx; + pgs_dev_ = pgs; +#endif + return 0; +} + +/* Create data for generators that is used in different computations */ +int GENParamsRajaHiop::allocate(OPFLOW opflow) { + PS ps = opflow->ps; + PetscInt loc, gloc = 0, geni = 0; + PSGEN gen; + PSBUS bus; + PetscInt i, j; + PetscErrorCode ierr; + + PetscFunctionBegin; + + /* Get the number of active generators (STATUS ON) */ + ierr = PSGetNumActiveGenerators(ps, &ngenON, NULL); + CHKERRQ(ierr); + + /* Allocate arrays on the host */ + auto &resmgr = umpire::ResourceManager::getInstance(); + h_allocator_ = resmgr.getAllocator("HOST"); + + cost_alpha = paramAlloc(h_allocator_, ngenON); + cost_beta = paramAlloc(h_allocator_, ngenON); + cost_gamma = paramAlloc(h_allocator_, ngenON); + + pt = paramAlloc(h_allocator_, ngenON); + pb = paramAlloc(h_allocator_, ngenON); + qt = paramAlloc(h_allocator_, ngenON); + qb = paramAlloc(h_allocator_, ngenON); + isrenewable = paramAlloc(h_allocator_, ngenON); + + xidx = paramAlloc(h_allocator_, ngenON); + gidxbus = paramAlloc(h_allocator_, ngenON); + + eqjacspbus_idx = paramAlloc(h_allocator_, ngenON); + eqjacsqbus_idx = paramAlloc(h_allocator_, ngenON); + hesssp_idx = paramAlloc(h_allocator_, ngenON); + + if (opflow->has_gensetpoint) { + geqidxgen = paramAlloc(h_allocator_, ngenON); + gineqidxgen = paramAlloc(h_allocator_, ngenON); + gbineqidxgen = paramAlloc(h_allocator_, ngenON); + eqjacspgen_idx = paramAlloc(h_allocator_, ngenON); + ineqjacspgen_idx = paramAlloc(h_allocator_, ngenON); + pgs = paramAlloc(h_allocator_, ngenON); + } + + /* Insert data in genparams */ + for (i = 0; i < ps->nbus; i++) { + bus = &ps->bus[i]; + gloc = bus->starteqloc; + + for (j = 0; j < bus->ngen; j++) { + ierr = PSBUSGetGen(bus, j, &gen); + CHKERRQ(ierr); + if (!gen->status) + continue; + + loc = gen->startxpowloc; + + cost_alpha[geni] = gen->cost_alpha; + cost_beta[geni] = gen->cost_beta; + cost_gamma[geni] = gen->cost_gamma; + pt[geni] = gen->pt; + pb[geni] = gen->pb; + qt[geni] = gen->qt; + qb[geni] = gen->qb; + isrenewable[geni] = (int)gen->isrenewable; + if (opflow->has_gensetpoint) { + pgs[geni] = gen->pgs; + } + + xidx[geni] = opflow->idxn2sd_map[loc]; + gidxbus[geni] = gloc; + if (opflow->has_gensetpoint) { + geqidxgen[geni] = gen->starteqloc; + gineqidxgen[geni] = gen->startineqloc; + gbineqidxgen[geni] = opflow->nconeq + gen->startineqloc; + } + + geni++; + } + } + +#ifdef EXAGO_ENABLE_GPU + d_allocator_ = resmgr.getAllocator("DEVICE"); + // Allocate arrays on the device + cost_alpha_dev_ = paramAlloc(d_allocator_, ngenON); + cost_beta_dev_ = paramAlloc(d_allocator_, ngenON); + cost_gamma_dev_ = paramAlloc(d_allocator_, ngenON); + + pt_dev_ = paramAlloc(d_allocator_, ngenON); + pb_dev_ = paramAlloc(d_allocator_, ngenON); + qt_dev_ = paramAlloc(d_allocator_, ngenON); + qb_dev_ = paramAlloc(d_allocator_, ngenON); + isrenewable_dev_ = paramAlloc(d_allocator_, ngenON); + + xidx_dev_ = paramAlloc(d_allocator_, ngenON); + gidxbus_dev_ = paramAlloc(d_allocator_, ngenON); + + eqjacspbus_idx_dev_ = paramAlloc(d_allocator_, ngenON); + eqjacsqbus_idx_dev_ = paramAlloc(d_allocator_, ngenON); + hesssp_idx_dev_ = paramAlloc(d_allocator_, ngenON); + if (opflow->has_gensetpoint) { + geqidxgen_dev_ = paramAlloc(d_allocator_, ngenON); + gineqidxgen_dev_ = paramAlloc(d_allocator_, ngenON); + gbineqidxgen_dev_ = paramAlloc(d_allocator_, ngenON); + eqjacspgen_idx_dev_ = paramAlloc(d_allocator_, ngenON); + ineqjacspgen_idx_dev_ = paramAlloc(d_allocator_, ngenON); + pgs_dev_ = paramAlloc(d_allocator_, ngenON); + } +#endif + return 0; +} + +void PbpolModelRajaHiop::destroy(OPFLOW opflow) { + loadparams.destroy(opflow); + lineparams.destroy(opflow); + busparams.destroy(opflow); + genparams.destroy(opflow); + +#ifdef EXAGO_ENABLE_GPU + if (i_jaceq != NULL) { + // These arrays get allocated only with sparse GPU model. For other models, + // they are not allocated. Hence, this business of checking if the array is + // NULL + + auto &resmgr = umpire::ResourceManager::getInstance(); + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + + h_allocator_.deallocate(i_jaceq); + h_allocator_.deallocate(j_jaceq); + h_allocator_.deallocate(val_jaceq); + + h_allocator_.deallocate(i_hess); + h_allocator_.deallocate(j_hess); + h_allocator_.deallocate(val_hess); + + if (opflow->nconineq) { + h_allocator_.deallocate(i_jacineq); + h_allocator_.deallocate(j_jacineq); + h_allocator_.deallocate(val_jacineq); + } + } +#endif +} diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.h b/src/opflow/model/power_bal_hiop/paramsrajahiop.h new file mode 100644 index 00000000..91f1fb9f --- /dev/null +++ b/src/opflow/model/power_bal_hiop/paramsrajahiop.h @@ -0,0 +1,260 @@ + +#include + +#pragma once + +#include +#include +#include + +/// Class storing bus parameters +struct BUSParamsRajaHiop { + // Host data + int nbus; /* Number of buses */ + int *isref; /* isref[i] = 1 if bus is reference bus */ + int *isisolated; /* isisolated[i] = 1 if bus is isolated bus */ + int *ispvpq; /* For all other buses */ + double *powerimbalance_penalty; /* Penalty for power imbalance */ + double *vmin; /* min. voltage magnitude limit */ + double *vmax; /* max. voltage magnitude limit */ + double *va; /* bus angle (from file only used in bounds) */ + double *vm; /* bus voltage magnitude (from file only used in bounds) */ + double *gl; /* bus shunt (conductance) */ + double *bl; /* bus shunt (suspectance) */ + int *xidx; /* starting locations for bus variables in X vector */ + int *xidxpimb; /* starting locations for power imbalance bus variables in X + vector */ + int *gidx; /* starting locations for bus balance equations in constraint + vector */ + int *jacsp_idx; /* Location number in the sparse Jacobian for Pimb */ + int *jacsq_idx; /* Location number in the sparse Jacobian for Qimb */ + int *hesssp_idx; /* Location number in the Hessian */ + + // Device data + int *isref_dev_; /* isref[i] = 1 if bus is reference bus */ + int *isisolated_dev_; /* isisolated[i] = 1 if bus is isolated bus */ + int *ispvpq_dev_; /* For all other buses */ + double *powerimbalance_penalty_dev_; /* Penalty for power imbalance */ + double *vmin_dev_; /* min. voltage magnitude limit */ + double *vmax_dev_; /* max. voltage magnitude limit */ + double *va_dev_; /* bus angle (from file only used in bounds) */ + double *vm_dev_; /* bus voltage magnitude (from file only used in bounds) */ + double *gl_dev_; /* bus shunt (conductance) */ + double *bl_dev_; /* bus shunt (suspectance) */ + int *xidx_dev_; /* starting locations for bus variables in X vector */ + int *xidxpimb_dev_; /* starting locations for power imbalance bus variables in + X vector */ + int *gidx_dev_; /* starting locations for bus balance equations in constraint + vector */ + int *jacsp_idx_dev_; /* Location number in the sparse Jacobian for Pimb */ + int *jacsq_idx_dev_; /* Location number in the sparse Jacobian for Qimb */ + int *hesssp_idx_dev_; /* Location number in the Hessian */ + + int allocate(OPFLOW); + int destroy(OPFLOW); + int copy(OPFLOW); + +private: + // Umpire memory allocators + umpire::Allocator h_allocator_; + umpire::Allocator d_allocator_; +}; + +struct GENParamsRajaHiop { +public: + // Host data + int ngenON; /* Number of generators with STATUS ON */ + double *cost_alpha; /* generator cost coefficients */ + double *cost_beta; /* generator cost coefficients */ + double *cost_gamma; /* generator cost coefficients */ + double *pt; /* min. active power gen. limits */ + double *pb; /* max. active power gen. limits */ + double *qt; /* min. reactive power gen. limits */ + double *qb; /* max. reactive power gen. limits */ + double *pgs; /* real power output setpoint */ + int *isrenewable; /* Is renewable generator? */ + + int *xidx; /* starting locations in X vector */ + int * + gidxbus; /* starting locations in constraint vector for bus constraints */ + int *geqidxgen; /* starting locations in equality constraint vector for gen + constraints */ + int *gineqidxgen; /* starting locations in inequality constraint vector for + gen constraints */ + int *gbineqidxgen; /* Starting location to insert contribution to inequality + constraint bound */ + + /* The following members are only used with HIOP */ + int *eqjacspbus_idx; /* Location number in the bus equality constraints sparse + Jacobian for Pg */ + int *eqjacsqbus_idx; /* Location number in the bus equality constraints sparse + Jacobian for Qg */ + int *eqjacspgen_idx; /* Location number in the gen equality constraints sparse + Jacobian for Pg */ + int *ineqjacspgen_idx; /* Location number in the bus equality constraints + sparse Jacobian for Pg */ + int *hesssp_idx; /* Location number in the Hessian */ + + // Device data + double *cost_alpha_dev_; /* generator cost coefficients */ + double *cost_beta_dev_; /* generator cost coefficients */ + double *cost_gamma_dev_; /* generator cost coefficients */ + double *pt_dev_; /* min. active power gen. limits */ + double *pb_dev_; /* max. active power gen. limits */ + double *qt_dev_; /* min. reactive power gen. limits */ + double *qb_dev_; /* max. reactive power gen. limits */ + double *pgs_dev_; /* real power output setpoint */ + int *isrenewable_dev_; /* Is renewable generator? */ + + int *xidx_dev_; /* starting locations in X vector */ + int *gidxbus_dev_; /* starting locations in constraint vector for bus + constraints */ + int *geqidxgen_dev_; /* starting locations in equality constraint vector for + gen constraints */ + int *gineqidxgen_dev_; /* starting locations in inequality constraint vector + for gen constraints */ + int *gbineqidxgen_dev_; /* Starting location to insert contribution to + inequality constraint bound */ + + /* The following members are only used with HIOP */ + int *eqjacspbus_idx_dev_; /* Location number in the bus equality constraints + sparse Jacobian for Pg */ + int *eqjacsqbus_idx_dev_; /* Location number in the bus equality constraints + sparse Jacobian for Qg */ + int *eqjacspgen_idx_dev_; /* Location number in the gen equality constraints + sparse Jacobian for Pg */ + int *ineqjacspgen_idx_dev_; /* Location number in the bus equality constraints + sparse Jacobian for Pg */ + int *hesssp_idx_dev_; /* Location number in the Hessian */ + + int allocate(OPFLOW); + int destroy(OPFLOW); + int copy(OPFLOW); + +private: + // Umpire memory allocators + umpire::Allocator h_allocator_; + umpire::Allocator d_allocator_; +}; + +struct LOADParamsRajaHiop { + // Host data + int nload; /* Number of loads */ + double *pl; /* active power demand */ + double *ql; /* reactive power demand */ + double *loadloss_penalty; /* Penalty for load loss */ + int *xidx; /* starting location in X vector */ + int *gidx; /* starting location in constraint vector */ + + /* The following members are only used with HIOP */ + int *jacsp_idx; /* Location number in the sparse Jacobian for delPload */ + int *jacsq_idx; /* Location number in the sparse Jacobian for delQload */ + int *hesssp_idx; /* Location number in the Hessian */ + + double *pl_dev_; /* active power demand */ + double *ql_dev_; /* reactive power demand */ + double *loadloss_penalty_dev_; /* Penalty for load loss */ + int *xidx_dev_; /* starting location in X vector */ + int *gidx_dev_; /* starting location in constraint vector */ + + int *jacsp_idx_dev_; /* Location number in the sparse Jacobian for delPload */ + int *jacsq_idx_dev_; /* Location number in the sparse Jacobian for delQload */ + int *hesssp_idx_dev_; /* Location number in the Hessian */ + + int allocate(OPFLOW); + int destroy(OPFLOW); + int copy(OPFLOW); + +private: + // Umpire memory allocators + umpire::Allocator h_allocator_; + umpire::Allocator d_allocator_; +}; + +struct LINEParamsRajaHiop { + // Host data + int nlineON; /* Number of active lines (STATUS = 1) */ + int nlinelim; /* Active lines + limits */ + double *Gff; /* From side self conductance */ + double *Bff; /* From side self susceptance */ + double *Gft; /* From-to transfer conductance */ + double *Bft; /* From-to transfer susceptance */ + double *Gtf; /* To-from transfer conductance */ + double *Btf; /* To-from transfer susceptance */ + double *Gtt; /* To side self conductance */ + double *Btt; /* To side self susceptance */ + double *rateA; /* Line MVA rating A (normal operation) */ + int *xidxf; /* Starting locatin of from bus voltage variables */ + int *xidxt; /* Starting location of to bus voltage variables */ + int *geqidxf; /* Starting location of from side to insert equality constraint + contribution in constraints vector */ + int *geqidxt; /* Starting location of to side to insert equality constraint + contribution in constraints vector */ + int *gineqidx; /* Starting location to insert contribution to inequality + constraint */ + int *gbineqidx; /* Starting location to insert contribution to inequality + constraint bound */ + int *linelimidx; /* Indices for subset of lines that have finite limits */ + + // Device data + double *Gff_dev_; /* From side self conductance */ + double *Bff_dev_; /* From side self susceptance */ + double *Gft_dev_; /* From-to transfer conductance */ + double *Bft_dev_; /* From-to transfer susceptance */ + double *Gtf_dev_; /* To-from transfer conductance */ + double *Btf_dev_; /* To-from transfer susceptance */ + double *Gtt_dev_; /* To side self conductance */ + double *Btt_dev_; /* To side self susceptance */ + double *rateA_dev_; /* Line MVA rating A (normal operation) */ + int *xidxf_dev_; /* Starting locatin of from bus voltage variables */ + int *xidxt_dev_; /* Starting location of to bus voltage variables */ + int *geqidxf_dev_; /* Starting location of from side to insert equality + constraint contribution in constraints vector */ + int *geqidxt_dev_; /* Starting location of to side to insert equality + constraint contribution in constraints vector */ + int *gineqidx_dev_; /* Starting location to insert contribution to inequality + constraint */ + int *gbineqidx_dev_; /* Starting location to insert contribution to inequality + constraint bound */ + int * + linelimidx_dev_; /* Indices for subset of lines that have finite limits */ + + int allocate(OPFLOW); + int destroy(OPFLOW); + int copy(OPFLOW); + +private: + // Umpire memory allocators + umpire::Allocator h_allocator_; + umpire::Allocator d_allocator_; +}; + +struct _p_FormPBPOLRAJAHIOP {}; +typedef struct _p_FormPBPOLRAJAHIOP *PBPOLRAJAHIOP; + +struct PbpolModelRajaHiop : public _p_FormPBPOLRAJAHIOP { + PbpolModelRajaHiop(void) { + i_jaceq = j_jaceq = i_jacineq = j_jacineq = NULL; + i_hess = j_hess = NULL; + val_jaceq = val_jacineq = val_hess = NULL; + } + + void destroy(OPFLOW opflow); + + ~PbpolModelRajaHiop() {} + + GENParamsRajaHiop genparams; + LOADParamsRajaHiop loadparams; + LINEParamsRajaHiop lineparams; + BUSParamsRajaHiop busparams; + + // Arrays to store Jacobian and Hessian indices and entries on CPU (used with + // GPU sparse model) + int *i_jaceq, + *j_jaceq; // Row and column indices for equality constrained Jacobian + int *i_jacineq, + *j_jacineq; // Row and column indices for inequality constrained Jacobain + int *i_hess, *j_hess; // Row and column indices for hessian + double *val_jaceq, *val_jacineq, + *val_hess; // values for equality, inequality jacobians and hessian +}; diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiop.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiop.cpp index 9226157a..0d7930bc 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiop.cpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiop.cpp @@ -509,7 +509,7 @@ extern PetscErrorCode OPFLOWModelSetNumConstraints_PBPOL(OPFLOW, PetscInt *, PetscErrorCode OPFLOWModelCreate_PBPOLRAJAHIOP(OPFLOW opflow) { PetscFunctionBegin; - PbpolModelRajaHiop *pbpol = new PbpolModelRajaHiop(opflow); + PbpolModelRajaHiop *pbpol = new PbpolModelRajaHiop(); opflow->model = pbpol; diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiop.h b/src/opflow/model/power_bal_hiop/pbpolrajahiop.h index 4a5bfac6..3092c548 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiop.h +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiop.h @@ -6,9 +6,7 @@ #define _PBPOLRAJAHIOP_H #include - -struct _p_FormPBPOLRAJAHIOP {}; -typedef struct _p_FormPBPOLRAJAHIOP *PBPOLRAJAHIOP; +#include "paramsrajahiop.h" extern PetscErrorCode OPFLOWSetVariableBoundsArray_PBPOLRAJAHIOP(OPFLOW, double *, double *); diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.cpp index d7e213df..0aa506bb 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.cpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.cpp @@ -17,767 +17,6 @@ /* No Load loss or power imbalance variables considered yet */ /********************************************/ -/* Functions to create and destroy data arrays for different - component classes -*/ -int BUSParamsRajaHiop::destroy(OPFLOW opflow) { - h_allocator_.deallocate(isref); - h_allocator_.deallocate(isisolated); - h_allocator_.deallocate(ispvpq); - h_allocator_.deallocate(vmin); - h_allocator_.deallocate(vmax); - h_allocator_.deallocate(va); - h_allocator_.deallocate(vm); - h_allocator_.deallocate(gl); - h_allocator_.deallocate(bl); - h_allocator_.deallocate(xidx); - h_allocator_.deallocate(gidx); - if (opflow->include_powerimbalance_variables) { - h_allocator_.deallocate(xidxpimb); - h_allocator_.deallocate(powerimbalance_penalty); - h_allocator_.deallocate(jacsp_idx); - h_allocator_.deallocate(jacsq_idx); - } - -#ifdef EXAGO_ENABLE_GPU - d_allocator_.deallocate(isref_dev_); - d_allocator_.deallocate(isisolated_dev_); - d_allocator_.deallocate(ispvpq_dev_); - d_allocator_.deallocate(vmin_dev_); - d_allocator_.deallocate(vmax_dev_); - d_allocator_.deallocate(va_dev_); - d_allocator_.deallocate(vm_dev_); - d_allocator_.deallocate(gl_dev_); - d_allocator_.deallocate(bl_dev_); - d_allocator_.deallocate(xidx_dev_); - d_allocator_.deallocate(gidx_dev_); - if (opflow->include_powerimbalance_variables) { - d_allocator_.deallocate(xidxpimb_dev_); - d_allocator_.deallocate(powerimbalance_penalty_dev_); - d_allocator_.deallocate(jacsp_idx_dev_); - d_allocator_.deallocate(jacsq_idx_dev_); - } -#endif - - return 0; -} - -/* Copy data from host to device */ -int BUSParamsRajaHiop::copy(OPFLOW opflow) { - /* Allocate the arrays */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - -#ifdef EXAGO_ENABLE_GPU - d_allocator_ = resmgr.getAllocator("DEVICE"); - - // Copy host data from the host to the device - resmgr.copy(isref_dev_, isref); - resmgr.copy(isisolated_dev_, isisolated); - resmgr.copy(ispvpq_dev_, ispvpq); - - resmgr.copy(vmin_dev_, vmin); - resmgr.copy(vmax_dev_, vmax); - resmgr.copy(va_dev_, va); - resmgr.copy(vm_dev_, vm); - resmgr.copy(gl_dev_, gl); - resmgr.copy(bl_dev_, bl); - - resmgr.copy(xidx_dev_, xidx); - resmgr.copy(gidx_dev_, gidx); - if (opflow->include_powerimbalance_variables) { - resmgr.copy(xidxpimb_dev_, xidxpimb); - resmgr.copy(jacsp_idx_dev_, jacsp_idx); - resmgr.copy(jacsq_idx_dev_, jacsq_idx); - resmgr.copy(powerimbalance_penalty_dev_, powerimbalance_penalty); - } -#else - isref_dev_ = isref; - isisolated_dev_ = isisolated; - ispvpq_dev_ = ispvpq; - vmin_dev_ = vmin; - vmax_dev_ = vmax; - va_dev_ = va; - vm_dev_ = vm; - gl_dev_ = gl; - bl_dev_ = bl; - xidx_dev_ = xidx; - xidxpimb_dev_ = xidxpimb; - gidx_dev_ = gidx; - jacsp_idx_dev_ = jacsp_idx; - jacsq_idx_dev_ = jacsq_idx; - powerimbalance_penalty_dev_ = powerimbalance_penalty; -#endif - return 0; -} - -/* Create data for buses that is used in different computations */ -int BUSParamsRajaHiop::allocate(OPFLOW opflow) { - PS ps = opflow->ps; - PetscInt loc; - PSBUS bus; - - nbus = ps->nbus; - /* Allocate the arrays */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - - // Allocate data on the host - isref = paramAlloc(h_allocator_, nbus); - isisolated = paramAlloc(h_allocator_, nbus); - ispvpq = paramAlloc(h_allocator_, nbus); - - vmin = paramAlloc(h_allocator_, nbus); - vmax = paramAlloc(h_allocator_, nbus); - va = paramAlloc(h_allocator_, nbus); - vm = paramAlloc(h_allocator_, nbus); - gl = paramAlloc(h_allocator_, nbus); - bl = paramAlloc(h_allocator_, nbus); - - xidx = paramAlloc(h_allocator_, nbus); - gidx = paramAlloc(h_allocator_, nbus); - - if (opflow->include_powerimbalance_variables) { - xidxpimb = paramAlloc(h_allocator_, nbus); - powerimbalance_penalty = paramAlloc(h_allocator_, nbus); - jacsp_idx = paramAlloc(h_allocator_, nbus); - jacsq_idx = paramAlloc(h_allocator_, nbus); - } - - /* Memzero arrays */ - resmgr.memset(isref, 0, nbus * sizeof(int)); - resmgr.memset(ispvpq, 0, nbus * sizeof(int)); - resmgr.memset(isisolated, 0, nbus * sizeof(int)); - - for (int i = 0; i < nbus; i++) { - bus = &ps->bus[i]; - loc = bus->startxVloc; - - xidx[i] = opflow->idxn2sd_map[loc]; - gidx[i] = bus->starteqloc; - - if (bus->ide == REF_BUS) - isref[i] = 1; - else if (bus->ide == ISOLATED_BUS) - isisolated[i] = 1; - else - ispvpq[i] = 1; - - if (opflow->genbusvoltagetype == FIXED_AT_SETPOINT) { - if (bus->ide == REF_BUS || bus->ide == PV_BUS) { - /* Hold voltage at reference and PV buses */ - vmin[i] = bus->vm; - vmax[i] = bus->vm; - } else { - vmin[i] = bus->Vmin; - vmax[i] = bus->Vmax; - } - } else { - vmin[i] = bus->Vmin; - vmax[i] = bus->Vmax; - } - vm[i] = bus->vm; - va[i] = bus->va; - gl[i] = bus->gl; - bl[i] = bus->bl; - - if (opflow->include_powerimbalance_variables) { - loc = bus->startxpimbloc; - xidxpimb[i] = opflow->idxn2sd_map[loc]; - powerimbalance_penalty[i] = opflow->powerimbalance_penalty; - } - } - -#ifdef EXAGO_ENABLE_GPU - d_allocator_ = resmgr.getAllocator("DEVICE"); - // Allocate data on the device - isref_dev_ = paramAlloc(d_allocator_, nbus); - isisolated_dev_ = paramAlloc(d_allocator_, nbus); - ispvpq_dev_ = paramAlloc(d_allocator_, nbus); - - vmin_dev_ = paramAlloc(d_allocator_, nbus); - vmax_dev_ = paramAlloc(d_allocator_, nbus); - va_dev_ = paramAlloc(d_allocator_, nbus); - vm_dev_ = paramAlloc(d_allocator_, nbus); - gl_dev_ = paramAlloc(d_allocator_, nbus); - bl_dev_ = paramAlloc(d_allocator_, nbus); - - xidx_dev_ = paramAlloc(d_allocator_, nbus); - gidx_dev_ = paramAlloc(d_allocator_, nbus); - - if (opflow->include_powerimbalance_variables) { - xidxpimb_dev_ = paramAlloc(d_allocator_, nbus); - powerimbalance_penalty_dev_ = paramAlloc(d_allocator_, nbus); - jacsp_idx_dev_ = paramAlloc(d_allocator_, nbus); - jacsq_idx_dev_ = paramAlloc(d_allocator_, nbus); - } -#endif - return 0; -} - -int LINEParamsRajaHiop::copy(OPFLOW opflow) { - /* Allocate the arrays */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - -#ifdef EXAGO_ENABLE_GPU - // Copy data from the host to the device - resmgr.copy(Gff_dev_, Gff); - resmgr.copy(Bff_dev_, Bff); - resmgr.copy(Gft_dev_, Gft); - resmgr.copy(Bft_dev_, Bft); - resmgr.copy(Gtf_dev_, Gtf); - resmgr.copy(Btf_dev_, Btf); - resmgr.copy(Gtt_dev_, Gtt); - resmgr.copy(Btt_dev_, Btt); - resmgr.copy(rateA_dev_, rateA); - - resmgr.copy(xidxf_dev_, xidxf); - resmgr.copy(xidxt_dev_, xidxt); - - resmgr.copy(geqidxf_dev_, geqidxf); - resmgr.copy(geqidxt_dev_, geqidxt); - - if (opflow->nlinesmon) { - resmgr.copy(gineqidx_dev_, gineqidx); - resmgr.copy(gbineqidx_dev_, gbineqidx); - resmgr.copy(linelimidx_dev_, linelimidx); - } -#else - Gff_dev_ = Gff; - Bff_dev_ = Bff; - Gft_dev_ = Gft; - Bft_dev_ = Bft; - Gtf_dev_ = Gtf; - Btf_dev_ = Btf; - Gtt_dev_ = Gtt; - Btt_dev_ = Btt; - rateA_dev_ = rateA; - xidxf_dev_ = xidxf; - xidxt_dev_ = xidxt; - geqidxf_dev_ = geqidxf; - geqidxt_dev_ = geqidxt; - if (opflow->nlinesmon) { - gineqidx_dev_ = gineqidx; - gbineqidx_dev_ = gbineqidx; - linelimidx_dev_ = linelimidx; - } -#endif - return 0; -} - -int LINEParamsRajaHiop::destroy(OPFLOW opflow) { - // Destroy parameter arrays on the host - h_allocator_.deallocate(Gff); - h_allocator_.deallocate(Bff); - h_allocator_.deallocate(Gft); - h_allocator_.deallocate(Bft); - h_allocator_.deallocate(Gtf); - h_allocator_.deallocate(Btf); - h_allocator_.deallocate(Gtt); - h_allocator_.deallocate(Btt); - h_allocator_.deallocate(rateA); - - h_allocator_.deallocate(xidxf); - h_allocator_.deallocate(xidxt); - - h_allocator_.deallocate(geqidxf); - h_allocator_.deallocate(geqidxt); - - if (opflow->nlinesmon) { - h_allocator_.deallocate(gineqidx); - h_allocator_.deallocate(gbineqidx); - h_allocator_.deallocate(linelimidx); - } - -#ifdef EXAGO_ENABLE_GPU - // Destroy parameter arrays on the device - d_allocator_.deallocate(Gff_dev_); - d_allocator_.deallocate(Bff_dev_); - d_allocator_.deallocate(Gft_dev_); - d_allocator_.deallocate(Bft_dev_); - d_allocator_.deallocate(Gtf_dev_); - d_allocator_.deallocate(Btf_dev_); - d_allocator_.deallocate(Gtt_dev_); - d_allocator_.deallocate(Btt_dev_); - d_allocator_.deallocate(rateA_dev_); - - d_allocator_.deallocate(xidxf_dev_); - d_allocator_.deallocate(xidxt_dev_); - - d_allocator_.deallocate(geqidxf_dev_); - d_allocator_.deallocate(geqidxt_dev_); - - if (opflow->nlinesmon) { - d_allocator_.deallocate(gineqidx_dev_); - d_allocator_.deallocate(gbineqidx_dev_); - d_allocator_.deallocate(linelimidx_dev_); - } -#endif - - return 0; -} - -/* Create data for lines that is used in different computations */ -int LINEParamsRajaHiop::allocate(OPFLOW opflow) { - PS ps = opflow->ps; - PetscInt linei = 0; - PSLINE line; - PetscInt i; - PetscErrorCode ierr; - const PSBUS *connbuses; - PSBUS busf, bust; - - PetscFunctionBegin; - ierr = PSGetNumActiveLines(ps, &nlineON, NULL); - CHKERRQ(ierr); - - nlinelim = opflow->nlinesmon; - - /* Allocate data arrays */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - - // Allocate data on the host - Gff = paramAlloc(h_allocator_, nlineON); - Bff = paramAlloc(h_allocator_, nlineON); - Gft = paramAlloc(h_allocator_, nlineON); - Bft = paramAlloc(h_allocator_, nlineON); - Gtf = paramAlloc(h_allocator_, nlineON); - Btf = paramAlloc(h_allocator_, nlineON); - Gtt = paramAlloc(h_allocator_, nlineON); - Btt = paramAlloc(h_allocator_, nlineON); - rateA = paramAlloc(h_allocator_, nlineON); - - xidxf = paramAlloc(h_allocator_, nlineON); - xidxt = paramAlloc(h_allocator_, nlineON); - - geqidxf = paramAlloc(h_allocator_, nlineON); - geqidxt = paramAlloc(h_allocator_, nlineON); - - if (opflow->nlinesmon) { - linelimidx = paramAlloc(h_allocator_, nlinelim); - gineqidx = paramAlloc(h_allocator_, nlinelim); - gbineqidx = paramAlloc(h_allocator_, nlinelim); - } - - PetscInt j = 0; - /* Populate arrays */ - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - - if (!line->status) - continue; - - Gff[linei] = line->yff[0]; - Bff[linei] = line->yff[1]; - Gft[linei] = line->yft[0]; - Bft[linei] = line->yft[1]; - Gtf[linei] = line->ytf[0]; - Btf[linei] = line->ytf[1]; - Gtt[linei] = line->ytt[0]; - Btt[linei] = line->ytt[1]; - rateA[linei] = line->rateA; - - ierr = PSLINEGetConnectedBuses(line, &connbuses); - CHKERRQ(ierr); - busf = connbuses[0]; - bust = connbuses[1]; - - int xidxfi, xidxti; - xidxfi = busf->startxVloc; - xidxti = bust->startxVloc; - - xidxf[linei] = opflow->idxn2sd_map[xidxfi]; - xidxt[linei] = opflow->idxn2sd_map[xidxti]; - - /* - Each bus has two equality (balance) constraints, hence the use of - coefficient 2 to map the location of the equality constraint for the bus - */ - geqidxf[linei] = busf->starteqloc; - geqidxt[linei] = bust->starteqloc; - - if (j < opflow->nlinesmon && opflow->linesmon[j] == i) { - gbineqidx[j] = opflow->nconeq + line->startineqloc; - gineqidx[j] = line->startineqloc; - linelimidx[j] = linei; - j++; - } - - linei++; - } - -#ifdef EXAGO_ENABLE_GPU - d_allocator_ = resmgr.getAllocator("DEVICE"); - // Allocate data on the device - Gff_dev_ = paramAlloc(d_allocator_, nlineON); - Bff_dev_ = paramAlloc(d_allocator_, nlineON); - Gft_dev_ = paramAlloc(d_allocator_, nlineON); - Bft_dev_ = paramAlloc(d_allocator_, nlineON); - Gtf_dev_ = paramAlloc(d_allocator_, nlineON); - Btf_dev_ = paramAlloc(d_allocator_, nlineON); - Gtt_dev_ = paramAlloc(d_allocator_, nlineON); - Btt_dev_ = paramAlloc(d_allocator_, nlineON); - rateA_dev_ = paramAlloc(d_allocator_, nlineON); - - xidxf_dev_ = paramAlloc(d_allocator_, nlineON); - xidxt_dev_ = paramAlloc(d_allocator_, nlineON); - - geqidxf_dev_ = paramAlloc(d_allocator_, nlineON); - geqidxt_dev_ = paramAlloc(d_allocator_, nlineON); - - if (opflow->nconineq) { - gineqidx_dev_ = paramAlloc(d_allocator_, nlinelim); - gbineqidx_dev_ = paramAlloc(d_allocator_, nlinelim); - linelimidx_dev_ = paramAlloc(d_allocator_, nlinelim); - } -#endif - return 0; -} - -int LOADParamsRajaHiop::copy(OPFLOW opflow) { - /* Allocate the arrays */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - -#ifdef EXAGO_ENABLE_GPU - // Copy data from host to device - resmgr.copy(pl_dev_, pl); - resmgr.copy(ql_dev_, ql); - resmgr.copy(xidx_dev_, xidx); - resmgr.copy(gidx_dev_, gidx); - if (opflow->include_loadloss_variables) { - resmgr.copy(jacsp_idx_dev_, jacsp_idx); - resmgr.copy(jacsq_idx_dev_, jacsq_idx); - resmgr.copy(hesssp_idx_dev_, hesssp_idx); - resmgr.copy(loadloss_penalty_dev_, loadloss_penalty); - } -#else - pl_dev_ = pl; - ql_dev_ = ql; - xidx_dev_ = xidx; - gidx_dev_ = gidx; - jacsp_idx_dev_ = jacsp_idx; - jacsq_idx_dev_ = jacsq_idx; - hesssp_idx_dev_ = hesssp_idx; - loadloss_penalty_dev_ = loadloss_penalty; -#endif - - return 0; -} - -int LOADParamsRajaHiop::destroy(OPFLOW opflow) { - h_allocator_.deallocate(pl); - h_allocator_.deallocate(ql); - h_allocator_.deallocate(xidx); - h_allocator_.deallocate(gidx); - if (opflow->include_loadloss_variables) { - h_allocator_.deallocate(loadloss_penalty); - h_allocator_.deallocate(jacsp_idx); - h_allocator_.deallocate(jacsq_idx); - h_allocator_.deallocate(hesssp_idx); - } -#ifdef EXAGO_ENABLE_GPU - d_allocator_.deallocate(pl_dev_); - d_allocator_.deallocate(ql_dev_); - d_allocator_.deallocate(xidx_dev_); - d_allocator_.deallocate(gidx_dev_); - if (opflow->include_loadloss_variables) { - d_allocator_.deallocate(loadloss_penalty_dev_); - d_allocator_.deallocate(jacsp_idx_dev_); - d_allocator_.deallocate(jacsq_idx_dev_); - d_allocator_.deallocate(hesssp_idx_dev_); - } -#endif - return 0; -} - -/* Create data for loads that is used in different computations */ -int LOADParamsRajaHiop::allocate(OPFLOW opflow) { - PS ps = opflow->ps; - PetscInt loc, loadi = 0; - PSLOAD load; - PSBUS bus; - PetscInt i, j; - PetscErrorCode ierr; - - /* Get the number of active generators (STATUS ON) */ - ierr = PSGetNumLoads(ps, &nload, NULL); - CHKERRQ(ierr); - - /* Allocate arrays */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - - // Allocate data on the host - pl = paramAlloc(h_allocator_, nload); - ql = paramAlloc(h_allocator_, nload); - xidx = paramAlloc(h_allocator_, nload); - gidx = paramAlloc(h_allocator_, nload); - if (opflow->include_loadloss_variables) { - loadloss_penalty = paramAlloc(h_allocator_, nload); - jacsp_idx = paramAlloc(h_allocator_, nload); - jacsq_idx = paramAlloc(h_allocator_, nload); - hesssp_idx = paramAlloc(h_allocator_, nload); - } - /* Insert data in loadparams */ - for (i = 0; i < ps->nbus; i++) { - bus = &ps->bus[i]; - ierr = PSBUSGetVariableLocation(bus, &loc); - CHKERRQ(ierr); - for (j = 0; j < bus->nload; j++) { - ierr = PSBUSGetLoad(bus, j, &load); - CHKERRQ(ierr); - pl[loadi] = load->pl; - ql[loadi] = load->ql; - if (opflow->include_loadloss_variables) { - loc = load->startxloadlossloc; - loadloss_penalty[loadi] = load->loss_cost; - xidx[loadi] = opflow->idxn2sd_map[loc]; - } - gidx[loadi] = bus->starteqloc; - loadi++; - } - } -#ifdef EXAGO_ENABLE_GPU - d_allocator_ = resmgr.getAllocator("DEVICE"); - // Allocate data on the device - pl_dev_ = paramAlloc(d_allocator_, nload); - ql_dev_ = paramAlloc(d_allocator_, nload); - xidx_dev_ = paramAlloc(d_allocator_, nload); - gidx_dev_ = paramAlloc(d_allocator_, nload); - if (opflow->include_loadloss_variables) { - loadloss_penalty_dev_ = paramAlloc(d_allocator_, nload); - jacsp_idx_dev_ = paramAlloc(d_allocator_, nload); - jacsq_idx_dev_ = paramAlloc(d_allocator_, nload); - hesssp_idx_dev_ = paramAlloc(d_allocator_, nload); - } - -#endif - - return (0); -} - -int GENParamsRajaHiop::destroy(OPFLOW opflow) { - // Free arrays on the host - h_allocator_.deallocate(cost_alpha); - h_allocator_.deallocate(cost_beta); - h_allocator_.deallocate(cost_gamma); - h_allocator_.deallocate(pt); - h_allocator_.deallocate(pb); - h_allocator_.deallocate(qt); - h_allocator_.deallocate(qb); - h_allocator_.deallocate(isrenewable); - h_allocator_.deallocate(xidx); - h_allocator_.deallocate(gidxbus); - h_allocator_.deallocate(eqjacspbus_idx); - h_allocator_.deallocate(eqjacsqbus_idx); - h_allocator_.deallocate(hesssp_idx); - if (opflow->has_gensetpoint) { - h_allocator_.deallocate(geqidxgen); - h_allocator_.deallocate(gineqidxgen); - h_allocator_.deallocate(gbineqidxgen); - h_allocator_.deallocate(pgs); - h_allocator_.deallocate(eqjacspgen_idx); - h_allocator_.deallocate(ineqjacspgen_idx); - } -#ifdef EXAGO_ENABLE_GPU - // Free arrays on the device - d_allocator_.deallocate(cost_alpha_dev_); - d_allocator_.deallocate(cost_beta_dev_); - d_allocator_.deallocate(cost_gamma_dev_); - d_allocator_.deallocate(pt_dev_); - d_allocator_.deallocate(pb_dev_); - d_allocator_.deallocate(qt_dev_); - d_allocator_.deallocate(qb_dev_); - d_allocator_.deallocate(isrenewable_dev_); - d_allocator_.deallocate(xidx_dev_); - d_allocator_.deallocate(gidxbus_dev_); - d_allocator_.deallocate(eqjacspbus_idx_dev_); - d_allocator_.deallocate(eqjacsqbus_idx_dev_); - d_allocator_.deallocate(hesssp_idx_dev_); - if (opflow->has_gensetpoint) { - d_allocator_.deallocate(geqidxgen_dev_); - d_allocator_.deallocate(gineqidxgen_dev_); - d_allocator_.deallocate(gbineqidxgen_dev_); - d_allocator_.deallocate(pgs_dev_); - d_allocator_.deallocate(eqjacspgen_idx_dev_); - d_allocator_.deallocate(ineqjacspgen_idx_dev_); - } - -#endif - return 0; -} - -int GENParamsRajaHiop::copy(OPFLOW opflow) { - /* Allocate arrays on the host */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - -#ifdef EXAGO_ENABLE_GPU - // Copy host data to the device - resmgr.copy(cost_alpha_dev_, cost_alpha); - resmgr.copy(cost_beta_dev_, cost_beta); - resmgr.copy(cost_gamma_dev_, cost_gamma); - - resmgr.copy(pt_dev_, pt); - resmgr.copy(pb_dev_, pb); - resmgr.copy(qt_dev_, qt); - resmgr.copy(qb_dev_, qb); - resmgr.copy(isrenewable_dev_, isrenewable); - - resmgr.copy(xidx_dev_, xidx); - resmgr.copy(gidxbus_dev_, gidxbus); - - resmgr.copy(eqjacspbus_idx_dev_, eqjacspbus_idx); - resmgr.copy(eqjacsqbus_idx_dev_, eqjacsqbus_idx); - resmgr.copy(hesssp_idx_dev_, hesssp_idx); - if (opflow->has_gensetpoint) { - resmgr.copy(geqidxgen_dev_, geqidxgen); - resmgr.copy(gineqidxgen_dev_, gineqidxgen); - resmgr.copy(gbineqidxgen_dev_, gbineqidxgen); - resmgr.copy(eqjacspgen_idx_dev_, eqjacspgen_idx); - resmgr.copy(ineqjacspgen_idx_dev_, ineqjacspgen_idx); - resmgr.copy(pgs_dev_, pgs); - } -#else - cost_alpha_dev_ = cost_alpha; - cost_beta_dev_ = cost_beta; - cost_gamma_dev_ = cost_gamma; - pt_dev_ = pt; - pb_dev_ = pb; - qt_dev_ = qt; - qb_dev_ = qb; - isrenewable_dev_ = isrenewable; - xidx_dev_ = xidx; - gidxbus_dev_ = gidxbus; - eqjacspbus_idx_dev_ = eqjacspbus_idx; - eqjacsqbus_idx_dev_ = eqjacsqbus_idx; - hesssp_idx_dev_ = hesssp_idx; - geqidxgen_dev_ = geqidxgen; - gineqidxgen_dev_ = gineqidxgen; - gbineqidxgen_dev_ = gbineqidxgen; - eqjacspgen_idx_dev_ = eqjacspgen_idx; - ineqjacspgen_idx_dev_ = ineqjacspgen_idx; - pgs_dev_ = pgs; -#endif - return 0; -} - -/* Create data for generators that is used in different computations */ -int GENParamsRajaHiop::allocate(OPFLOW opflow) { - PS ps = opflow->ps; - PetscInt loc, gloc = 0, geni = 0; - PSGEN gen; - PSBUS bus; - PetscInt i, j; - PetscErrorCode ierr; - - PetscFunctionBegin; - - /* Get the number of active generators (STATUS ON) */ - ierr = PSGetNumActiveGenerators(ps, &ngenON, NULL); - CHKERRQ(ierr); - - /* Allocate arrays on the host */ - auto &resmgr = umpire::ResourceManager::getInstance(); - h_allocator_ = resmgr.getAllocator("HOST"); - - cost_alpha = paramAlloc(h_allocator_, ngenON); - cost_beta = paramAlloc(h_allocator_, ngenON); - cost_gamma = paramAlloc(h_allocator_, ngenON); - - pt = paramAlloc(h_allocator_, ngenON); - pb = paramAlloc(h_allocator_, ngenON); - qt = paramAlloc(h_allocator_, ngenON); - qb = paramAlloc(h_allocator_, ngenON); - isrenewable = paramAlloc(h_allocator_, ngenON); - - xidx = paramAlloc(h_allocator_, ngenON); - gidxbus = paramAlloc(h_allocator_, ngenON); - - eqjacspbus_idx = paramAlloc(h_allocator_, ngenON); - eqjacsqbus_idx = paramAlloc(h_allocator_, ngenON); - hesssp_idx = paramAlloc(h_allocator_, ngenON); - - if (opflow->has_gensetpoint) { - geqidxgen = paramAlloc(h_allocator_, ngenON); - gineqidxgen = paramAlloc(h_allocator_, ngenON); - gbineqidxgen = paramAlloc(h_allocator_, ngenON); - eqjacspgen_idx = paramAlloc(h_allocator_, ngenON); - ineqjacspgen_idx = paramAlloc(h_allocator_, ngenON); - pgs = paramAlloc(h_allocator_, ngenON); - } - - /* Insert data in genparams */ - for (i = 0; i < ps->nbus; i++) { - bus = &ps->bus[i]; - gloc = bus->starteqloc; - - for (j = 0; j < bus->ngen; j++) { - ierr = PSBUSGetGen(bus, j, &gen); - CHKERRQ(ierr); - if (!gen->status) - continue; - - loc = gen->startxpowloc; - - cost_alpha[geni] = gen->cost_alpha; - cost_beta[geni] = gen->cost_beta; - cost_gamma[geni] = gen->cost_gamma; - pt[geni] = gen->pt; - pb[geni] = gen->pb; - qt[geni] = gen->qt; - qb[geni] = gen->qb; - isrenewable[geni] = (int)gen->isrenewable; - if (opflow->has_gensetpoint) { - pgs[geni] = gen->pgs; - } - - xidx[geni] = opflow->idxn2sd_map[loc]; - gidxbus[geni] = gloc; - if (opflow->has_gensetpoint) { - geqidxgen[geni] = gen->starteqloc; - gineqidxgen[geni] = gen->startineqloc; - gbineqidxgen[geni] = opflow->nconeq + gen->startineqloc; - } - - geni++; - } - } - -#ifdef EXAGO_ENABLE_GPU - d_allocator_ = resmgr.getAllocator("DEVICE"); - // Allocate arrays on the device - cost_alpha_dev_ = paramAlloc(d_allocator_, ngenON); - cost_beta_dev_ = paramAlloc(d_allocator_, ngenON); - cost_gamma_dev_ = paramAlloc(d_allocator_, ngenON); - - pt_dev_ = paramAlloc(d_allocator_, ngenON); - pb_dev_ = paramAlloc(d_allocator_, ngenON); - qt_dev_ = paramAlloc(d_allocator_, ngenON); - qb_dev_ = paramAlloc(d_allocator_, ngenON); - isrenewable_dev_ = paramAlloc(d_allocator_, ngenON); - - xidx_dev_ = paramAlloc(d_allocator_, ngenON); - gidxbus_dev_ = paramAlloc(d_allocator_, ngenON); - - eqjacspbus_idx_dev_ = paramAlloc(d_allocator_, ngenON); - eqjacsqbus_idx_dev_ = paramAlloc(d_allocator_, ngenON); - hesssp_idx_dev_ = paramAlloc(d_allocator_, ngenON); - if (opflow->has_gensetpoint) { - geqidxgen_dev_ = paramAlloc(d_allocator_, ngenON); - gineqidxgen_dev_ = paramAlloc(d_allocator_, ngenON); - gbineqidxgen_dev_ = paramAlloc(d_allocator_, ngenON); - eqjacspgen_idx_dev_ = paramAlloc(d_allocator_, ngenON); - ineqjacspgen_idx_dev_ = paramAlloc(d_allocator_, ngenON); - pgs_dev_ = paramAlloc(d_allocator_, ngenON); - } -#endif - return 0; -} - PetscErrorCode OPFLOWSetInitialGuessArray_PBPOLRAJAHIOP(OPFLOW opflow, double *x0_dev) { PetscErrorCode ierr; diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.h b/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.h index 5c36563c..682089a6 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.h +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopkernels.h @@ -4,11 +4,7 @@ #pragma once -#include -#include - #include "pbpolrajahiop.h" -#include /** @@ -28,248 +24,4 @@ void registerWith(T *ptr, SizeType N, umpire::ResourceManager &resmgr, resmgr.registerAllocation(ptr, record); } -template T *paramAlloc(umpire::Allocator &a, I N) { - return static_cast(a.allocate(N * sizeof(T))); -} - -/// Class storing bus parameters -struct BUSParamsRajaHiop { - // Host data - int nbus; /* Number of buses */ - int *isref; /* isref[i] = 1 if bus is reference bus */ - int *isisolated; /* isisolated[i] = 1 if bus is isolated bus */ - int *ispvpq; /* For all other buses */ - double *powerimbalance_penalty; /* Penalty for power imbalance */ - double *vmin; /* min. voltage magnitude limit */ - double *vmax; /* max. voltage magnitude limit */ - double *va; /* bus angle (from file only used in bounds) */ - double *vm; /* bus voltage magnitude (from file only used in bounds) */ - double *gl; /* bus shunt (conductance) */ - double *bl; /* bus shunt (suspectance) */ - int *xidx; /* starting locations for bus variables in X vector */ - int *xidxpimb; /* starting locations for power imbalance bus variables in X - vector */ - int *gidx; /* starting locations for bus balance equations in constraint - vector */ - int *jacsp_idx; /* Location number in the sparse Jacobian for Pimb */ - int *jacsq_idx; /* Location number in the sparse Jacobian for Qimb */ - int *hesssp_idx; /* Location number in the Hessian */ - - // Device data - int *isref_dev_; /* isref[i] = 1 if bus is reference bus */ - int *isisolated_dev_; /* isisolated[i] = 1 if bus is isolated bus */ - int *ispvpq_dev_; /* For all other buses */ - double *powerimbalance_penalty_dev_; /* Penalty for power imbalance */ - double *vmin_dev_; /* min. voltage magnitude limit */ - double *vmax_dev_; /* max. voltage magnitude limit */ - double *va_dev_; /* bus angle (from file only used in bounds) */ - double *vm_dev_; /* bus voltage magnitude (from file only used in bounds) */ - double *gl_dev_; /* bus shunt (conductance) */ - double *bl_dev_; /* bus shunt (suspectance) */ - int *xidx_dev_; /* starting locations for bus variables in X vector */ - int *xidxpimb_dev_; /* starting locations for power imbalance bus variables in - X vector */ - int *gidx_dev_; /* starting locations for bus balance equations in constraint - vector */ - int *jacsp_idx_dev_; /* Location number in the sparse Jacobian for Pimb */ - int *jacsq_idx_dev_; /* Location number in the sparse Jacobian for Qimb */ - int *hesssp_idx_dev_; /* Location number in the Hessian */ - - int allocate(OPFLOW); - int destroy(OPFLOW); - int copy(OPFLOW); - -private: - // Umpire memory allocators - umpire::Allocator h_allocator_; - umpire::Allocator d_allocator_; -}; - -struct GENParamsRajaHiop { -public: - // Host data - int ngenON; /* Number of generators with STATUS ON */ - double *cost_alpha; /* generator cost coefficients */ - double *cost_beta; /* generator cost coefficients */ - double *cost_gamma; /* generator cost coefficients */ - double *pt; /* min. active power gen. limits */ - double *pb; /* max. active power gen. limits */ - double *qt; /* min. reactive power gen. limits */ - double *qb; /* max. reactive power gen. limits */ - double *pgs; /* real power output setpoint */ - int *isrenewable; /* Is renewable generator? */ - - int *xidx; /* starting locations in X vector */ - int * - gidxbus; /* starting locations in constraint vector for bus constraints */ - int *geqidxgen; /* starting locations in equality constraint vector for gen - constraints */ - int *gineqidxgen; /* starting locations in inequality constraint vector for - gen constraints */ - int *gbineqidxgen; /* Starting location to insert contribution to inequality - constraint bound */ - - /* The following members are only used with HIOP */ - int *eqjacspbus_idx; /* Location number in the bus equality constraints sparse - Jacobian for Pg */ - int *eqjacsqbus_idx; /* Location number in the bus equality constraints sparse - Jacobian for Qg */ - int *eqjacspgen_idx; /* Location number in the gen equality constraints sparse - Jacobian for Pg */ - int *ineqjacspgen_idx; /* Location number in the bus equality constraints - sparse Jacobian for Pg */ - int *hesssp_idx; /* Location number in the Hessian */ - - // Device data - double *cost_alpha_dev_; /* generator cost coefficients */ - double *cost_beta_dev_; /* generator cost coefficients */ - double *cost_gamma_dev_; /* generator cost coefficients */ - double *pt_dev_; /* min. active power gen. limits */ - double *pb_dev_; /* max. active power gen. limits */ - double *qt_dev_; /* min. reactive power gen. limits */ - double *qb_dev_; /* max. reactive power gen. limits */ - double *pgs_dev_; /* real power output setpoint */ - int *isrenewable_dev_; /* Is renewable generator? */ - - int *xidx_dev_; /* starting locations in X vector */ - int *gidxbus_dev_; /* starting locations in constraint vector for bus - constraints */ - int *geqidxgen_dev_; /* starting locations in equality constraint vector for - gen constraints */ - int *gineqidxgen_dev_; /* starting locations in inequality constraint vector - for gen constraints */ - int *gbineqidxgen_dev_; /* Starting location to insert contribution to - inequality constraint bound */ - - /* The following members are only used with HIOP */ - int *eqjacspbus_idx_dev_; /* Location number in the bus equality constraints - sparse Jacobian for Pg */ - int *eqjacsqbus_idx_dev_; /* Location number in the bus equality constraints - sparse Jacobian for Qg */ - int *eqjacspgen_idx_dev_; /* Location number in the gen equality constraints - sparse Jacobian for Pg */ - int *ineqjacspgen_idx_dev_; /* Location number in the bus equality constraints - sparse Jacobian for Pg */ - int *hesssp_idx_dev_; /* Location number in the Hessian */ - - int allocate(OPFLOW); - int destroy(OPFLOW); - int copy(OPFLOW); - -private: - // Umpire memory allocators - umpire::Allocator h_allocator_; - umpire::Allocator d_allocator_; -}; - -struct LOADParamsRajaHiop { - // Host data - int nload; /* Number of loads */ - double *pl; /* active power demand */ - double *ql; /* reactive power demand */ - double *loadloss_penalty; /* Penalty for load loss */ - int *xidx; /* starting location in X vector */ - int *gidx; /* starting location in constraint vector */ - - /* The following members are only used with HIOP */ - int *jacsp_idx; /* Location number in the sparse Jacobian for delPload */ - int *jacsq_idx; /* Location number in the sparse Jacobian for delQload */ - int *hesssp_idx; /* Location number in the Hessian */ - - double *pl_dev_; /* active power demand */ - double *ql_dev_; /* reactive power demand */ - double *loadloss_penalty_dev_; /* Penalty for load loss */ - int *xidx_dev_; /* starting location in X vector */ - int *gidx_dev_; /* starting location in constraint vector */ - - int *jacsp_idx_dev_; /* Location number in the sparse Jacobian for delPload */ - int *jacsq_idx_dev_; /* Location number in the sparse Jacobian for delQload */ - int *hesssp_idx_dev_; /* Location number in the Hessian */ - - int allocate(OPFLOW); - int destroy(OPFLOW); - int copy(OPFLOW); - -private: - // Umpire memory allocators - umpire::Allocator h_allocator_; - umpire::Allocator d_allocator_; -}; - -struct LINEParamsRajaHiop { - // Host data - int nlineON; /* Number of active lines (STATUS = 1) */ - int nlinelim; /* Active lines + limits */ - double *Gff; /* From side self conductance */ - double *Bff; /* From side self susceptance */ - double *Gft; /* From-to transfer conductance */ - double *Bft; /* From-to transfer susceptance */ - double *Gtf; /* To-from transfer conductance */ - double *Btf; /* To-from transfer susceptance */ - double *Gtt; /* To side self conductance */ - double *Btt; /* To side self susceptance */ - double *rateA; /* Line MVA rating A (normal operation) */ - int *xidxf; /* Starting locatin of from bus voltage variables */ - int *xidxt; /* Starting location of to bus voltage variables */ - int *geqidxf; /* Starting location of from side to insert equality constraint - contribution in constraints vector */ - int *geqidxt; /* Starting location of to side to insert equality constraint - contribution in constraints vector */ - int *gineqidx; /* Starting location to insert contribution to inequality - constraint */ - int *gbineqidx; /* Starting location to insert contribution to inequality - constraint bound */ - int *linelimidx; /* Indices for subset of lines that have finite limits */ - - // Device data - double *Gff_dev_; /* From side self conductance */ - double *Bff_dev_; /* From side self susceptance */ - double *Gft_dev_; /* From-to transfer conductance */ - double *Bft_dev_; /* From-to transfer susceptance */ - double *Gtf_dev_; /* To-from transfer conductance */ - double *Btf_dev_; /* To-from transfer susceptance */ - double *Gtt_dev_; /* To side self conductance */ - double *Btt_dev_; /* To side self susceptance */ - double *rateA_dev_; /* Line MVA rating A (normal operation) */ - int *xidxf_dev_; /* Starting locatin of from bus voltage variables */ - int *xidxt_dev_; /* Starting location of to bus voltage variables */ - int *geqidxf_dev_; /* Starting location of from side to insert equality - constraint contribution in constraints vector */ - int *geqidxt_dev_; /* Starting location of to side to insert equality - constraint contribution in constraints vector */ - int *gineqidx_dev_; /* Starting location to insert contribution to inequality - constraint */ - int *gbineqidx_dev_; /* Starting location to insert contribution to inequality - constraint bound */ - int * - linelimidx_dev_; /* Indices for subset of lines that have finite limits */ - - int allocate(OPFLOW); - int destroy(OPFLOW); - int copy(OPFLOW); - -private: - // Umpire memory allocators - umpire::Allocator h_allocator_; - umpire::Allocator d_allocator_; -}; - -struct PbpolModelRajaHiop : public _p_FormPBPOLRAJAHIOP { - PbpolModelRajaHiop(OPFLOW opflow) { (void)opflow; } - - void destroy(OPFLOW opflow) { - loadparams.destroy(opflow); - lineparams.destroy(opflow); - busparams.destroy(opflow); - genparams.destroy(opflow); - } - - ~PbpolModelRajaHiop() {} - - GENParamsRajaHiop genparams; - LOADParamsRajaHiop loadparams; - LINEParamsRajaHiop lineparams; - BUSParamsRajaHiop busparams; -}; - #endif // EXAGO_ENABLE_HIOP diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp new file mode 100644 index 00000000..fa609378 --- /dev/null +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp @@ -0,0 +1,325 @@ +#include + +#if defined(EXAGO_ENABLE_RAJA) +#if defined(EXAGO_ENABLE_HIOP_SPARSE) + +#include +#include "pbpolrajahiopsparsekernels.hpp" + +/* Initialization is done on the host through this function. Copying over values + * to the device is done in OPFLOWSetInitialGuessArray_PBPOLRAJAHIOPSPARSE + */ +extern PetscErrorCode OPFLOWSetInitialGuess_PBPOL(OPFLOW, Vec, Vec); + +PetscErrorCode OPFLOWSetInitialGuess_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, Vec X, + Vec Lambda) { + PetscErrorCode ierr; + + PetscFunctionBegin; + + ierr = OPFLOWSetInitialGuess_PBPOL(opflow, X, Lambda); + CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +extern PetscErrorCode OPFLOWSetConstraintBounds_PBPOL(OPFLOW, Vec, Vec); + +/* The constraint bounds are also calculated on the host. + */ +PetscErrorCode OPFLOWSetConstraintBounds_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, + Vec Gl, Vec Gu) { + PetscErrorCode ierr; + + PetscFunctionBegin; + ierr = OPFLOWSetConstraintBounds_PBPOL(opflow, Gl, Gu); + CHKERRQ(ierr); + PetscFunctionReturn(0); +} + +extern PetscErrorCode OPFLOWSetVariableBounds_PBPOL(OPFLOW, Vec, Vec); + +/* The variable bounds are also calculated on the host. + */ +PetscErrorCode OPFLOWSetVariableBounds_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, + Vec Xl, Vec Xu) { + PetscErrorCode ierr; + + PetscFunctionBegin; + ierr = OPFLOWSetVariableBounds_PBPOL(opflow, Xl, Xu); + CHKERRQ(ierr); + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolutionToPS_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { + PetscErrorCode ierr; + PS ps = (PS)opflow->ps; + PetscInt i, k; + Vec X, Lambda; + PSBUS bus; + PSGEN gen; + PSLOAD load; + PSLINE line; + const PetscScalar *x, *lambda, *lambdae, *lambdai; + PetscInt loc, gloc = 0; + PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; + PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf; + PetscScalar Pf, Qf, Pt, Qt; + PSBUS busf, bust; + const PSBUS *connbuses; + PetscInt xlocf, xloct; + + PetscFunctionBegin; + + ierr = OPFLOWGetSolution(opflow, &X); + CHKERRQ(ierr); + ierr = OPFLOWGetConstraintMultipliers(opflow, &Lambda); + CHKERRQ(ierr); + + ierr = VecGetArrayRead(X, &x); + CHKERRQ(ierr); + ierr = VecGetArrayRead(Lambda, &lambda); + CHKERRQ(ierr); + lambdae = lambda; + if (opflow->Nconineq) { + lambdai = lambdae + opflow->nconeq; + } + + for (i = 0; i < ps->nbus; i++) { + bus = &ps->bus[i]; + + loc = bus->startxVloc; + + bus->va = x[loc]; + bus->vm = x[loc + 1]; + + gloc = bus->starteqloc; + bus->mult_pmis = lambdae[gloc]; + bus->mult_qmis = lambdae[gloc + 1]; + + if (opflow->include_powerimbalance_variables) { + loc = bus->startxpimbloc; + bus->pimb = x[loc]; + bus->qimb = x[loc + 1]; + } + + for (k = 0; k < bus->ngen; k++) { + ierr = PSBUSGetGen(bus, k, &gen); + CHKERRQ(ierr); + if (!gen->status) { + gen->pg = gen->qg = 0.0; + continue; + } + loc = gen->startxpowloc; + + gen->pg = x[loc]; + gen->qg = x[loc + 1]; + + if (opflow->has_gensetpoint) { + gloc += gen->nconeq; + } + } + + if (opflow->include_loadloss_variables) { + for (k = 0; k < bus->nload; k++) { + ierr = PSBUSGetLoad(bus, k, &load); + CHKERRQ(ierr); + loc = load->startxloadlossloc; + load->pl = load->pl - x[loc]; + load->ql = load->ql - x[loc + 1]; + } + } + } + + if (!opflow->ignore_lineflow_constraints) { + for (i = 0; i < ps->nline; i++) { + line = &ps->line[i]; + if (!line->status) { + line->mult_sf = line->mult_st = 0.0; + continue; + } + + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; + + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; + + xlocf = busf->startxVloc; + xloct = bust->startxVloc; + + thetaf = x[xlocf]; + Vmf = x[xlocf + 1]; + thetat = x[xloct]; + Vmt = x[xloct + 1]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; + + Pf = Gff * Vmf * Vmf + + Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + Qf = -Bff * Vmf * Vmf + + Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + + Pt = Gtt * Vmt * Vmt + + Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + Qt = -Btt * Vmt * Vmt + + Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + + line->pf = Pf; + line->qf = Qf; + line->pt = Pt; + line->qt = Qt; + line->sf = PetscSqrtScalar(Pf * Pf + Qf * Qf); + line->st = PetscSqrtScalar(Pt * Pt + Qt * Qt); + + if (line->rateA > 1e5) { + line->mult_sf = line->mult_st = 0.0; + } else { + gloc = line->startineqloc; + line->mult_sf = lambdai[gloc]; + line->mult_st = lambdai[gloc + 1]; + } + } + } + + ierr = VecRestoreArrayRead(X, &x); + CHKERRQ(ierr); + ierr = VecRestoreArrayRead(Lambda, &lambda); + CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +/* Reuse PBPOL model set up for obtaining locations */ +extern PetscErrorCode OPFLOWModelSetUp_PBPOL(OPFLOW); + +PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { + PetscErrorCode ierr; + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + + PetscFunctionBegin; + + ierr = OPFLOWModelSetUp_PBPOL(opflow); + CHKERRQ(ierr); + + ierr = pbpolrajahiopsparse->busparams.allocate(opflow); + ierr = pbpolrajahiopsparse->genparams.allocate(opflow); + ierr = pbpolrajahiopsparse->lineparams.allocate(opflow); + ierr = pbpolrajahiopsparse->loadparams.allocate(opflow); + + BUSParamsRajaHiop *busparams = &pbpolrajahiopsparse->busparams; + GENParamsRajaHiop *genparams = &pbpolrajahiopsparse->genparams; + LOADParamsRajaHiop *loadparams = &pbpolrajahiopsparse->loadparams; + LINEParamsRajaHiop *lineparams = &pbpolrajahiopsparse->lineparams; + + /* Need to compute the number of nonzeros in equality, inequality constraint + * Jacobians and Hessian */ + int nnz_eqjacsp = 0, nnz_ineqjacsp = 0, nnz_hesssp = 0; + opflow->nnz_eqjacsp = nnz_eqjacsp; + opflow->nnz_ineqjacsp = nnz_ineqjacsp; + opflow->nnz_hesssp = nnz_hesssp; + + ierr = busparams->copy(opflow); + ierr = genparams->copy(opflow); + ierr = lineparams->copy(opflow); + ierr = loadparams->copy(opflow); + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWModelDestroy_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + + PetscFunctionBegin; + pbpolrajahiopsparse->destroy(opflow); + delete pbpolrajahiopsparse; + pbpolrajahiopsparse = nullptr; + + PetscFunctionReturn(0); +} + +/* reuse numvariables and numconstraints functions from PBPOL model */ +extern PetscErrorCode OPFLOWModelSetNumVariables_PBPOL(OPFLOW, PetscInt *, + PetscInt *, PetscInt *); +extern PetscErrorCode OPFLOWModelSetNumConstraints_PBPOL(OPFLOW, PetscInt *, + PetscInt *, PetscInt *, + PetscInt *); +extern PetscErrorCode OPFLOWComputeEqualityConstraintJacobian_PBPOL(OPFLOW, Vec, + Mat); +extern PetscErrorCode OPFLOWComputeInequalityConstraintJacobian_PBPOL(OPFLOW, + Vec, Mat); +extern PetscErrorCode OPFLOWComputeHessian_PBPOL(OPFLOW, Vec, Vec, Vec, Mat); + +extern PetscErrorCode OPFLOWSolutionCallback_PBPOLRAJAHIOPSPARSE( + OPFLOW, const double *, const double *, const double *, const double *, + const double *, double); + +PetscErrorCode OPFLOWModelCreate_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { + + PetscFunctionBegin; + + PbpolModelRajaHiop *pbpolrajahiopsparse = new PbpolModelRajaHiop(); + + opflow->model = pbpolrajahiopsparse; + + /* PBPOLRAJAHIOPSPARSE models only support VARIABLE_WITHIN_BOUNDS + * opflow->genbusvoltagetype + */ + opflow->genbusvoltagetype = VARIABLE_WITHIN_BOUNDS; + + opflow->spdnordering = PETSC_FALSE; + + /* Inherit Ops */ + opflow->modelops.destroy = OPFLOWModelDestroy_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setnumvariables = OPFLOWModelSetNumVariables_PBPOL; + opflow->modelops.setnumconstraints = OPFLOWModelSetNumConstraints_PBPOL; + opflow->modelops.setvariablebounds = + OPFLOWSetVariableBounds_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setvariableboundsarray = + OPFLOWSetVariableBoundsArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setconstraintbounds = + OPFLOWSetConstraintBounds_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setconstraintboundsarray = + OPFLOWSetConstraintBoundsArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setinitialguess = OPFLOWSetInitialGuess_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setinitialguessarray = + OPFLOWSetInitialGuessArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computeequalityconstraintsarray = + OPFLOWComputeEqualityConstraintsArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computeinequalityconstraintsarray = + OPFLOWComputeInequalityConstraintsArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computeobjectivearray = + OPFLOWComputeObjectiveArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computegradientarray = + OPFLOWComputeGradientArray_PBPOLRAJAHIOPSPARSE; + opflow->modelops.solutiontops = OPFLOWSolutionToPS_PBPOLRAJAHIOPSPARSE; + opflow->modelops.setup = OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computeequalityconstraintjacobian = + OPFLOWComputeEqualityConstraintJacobian_PBPOL; + opflow->modelops.computesparseequalityconstraintjacobianhiop = + OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computeinequalityconstraintjacobian = + OPFLOWComputeInequalityConstraintJacobian_PBPOL; + opflow->modelops.computesparseinequalityconstraintjacobianhiop = + OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE; + opflow->modelops.computehessian = OPFLOWComputeHessian_PBPOL; + opflow->modelops.computesparsehessianhiop = + OPFLOWComputeSparseHessian_PBPOLRAJAHIOPSPARSE; + opflow->modelops.solutioncallbackhiop = + OPFLOWSolutionCallback_PBPOLRAJAHIOPSPARSE; + + PetscFunctionReturn(0); +} + +#endif +#endif diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.hpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.hpp new file mode 100644 index 00000000..ad549d69 --- /dev/null +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.hpp @@ -0,0 +1,51 @@ +#include + +#if defined(EXAGO_ENABLE_RAJA) +#if defined(EXAGO_ENABLE_HIOP_SPARSE) + +#ifndef _PBPOLRAJAHIOPSPARSE_H +#define _PBPOLRAJAHIOPSPARSE_H + +#include "pbpolrajahiop.h" +#include "paramsrajahiop.h" + +extern PetscErrorCode +OPFLOWSetVariableBoundsArray_PBPOLRAJAHIOPSPARSE(OPFLOW, double *, double *); +extern PetscErrorCode +OPFLOWSetConstraintBoundsArray_PBPOLRAJAHIOPSPARSE(OPFLOW, double *, double *); +extern PetscErrorCode OPFLOWSetInitialGuessArray_PBPOLRAJAHIOPSPARSE(OPFLOW, + double *); +extern PetscErrorCode OPFLOWComputeEqualityConstraintsArray_PBPOLRAJAHIOPSPARSE( + OPFLOW, const double *, double *); +extern PetscErrorCode +OPFLOWComputeInequalityConstraintsArray_PBPOLRAJAHIOPSPARSE(OPFLOW, + const double *, + double *); +extern PetscErrorCode +OPFLOWComputeObjectiveArray_PBPOLRAJAHIOPSPARSE(OPFLOW, const double *, + double *); +extern PetscErrorCode +OPFLOWComputeGradientArray_PBPOLRAJAHIOPSPARSE(OPFLOW, const double *, + double *); +extern PetscErrorCode +OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( + OPFLOW, const double *, int *, int *, double *); +extern PetscErrorCode +OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( + OPFLOW, const double *, int *, int *, double *); +extern PetscErrorCode OPFLOWComputeSparseHessian_PBPOLRAJAHIOPSPARSE( + OPFLOW, const double *, const double *, int *, int *, double *); +extern PetscErrorCode +OPFLOWComputeDenseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE(OPFLOW, + const double *, + double *); +extern PetscErrorCode +OPFLOWComputeDenseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( + OPFLOW, const double *, double *); +extern PetscErrorCode +OPFLOWComputeDenseHessian_PBPOLRAJAHIOPSPARSE(OPFLOW, const double *, + const double *, double *); + +#endif +#endif +#endif diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp new file mode 100644 index 00000000..adc10c8e --- /dev/null +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp @@ -0,0 +1,932 @@ + +#include + +#if defined(EXAGO_ENABLE_RAJA) +#if defined(EXAGO_ENABLE_HIOP_SPARSE) + +#include +#include + +#include +#include + +#include +#include +#include "pbpolrajahiopsparsekernels.hpp" +#include "pbpolrajahiopsparse.hpp" + +PetscErrorCode OPFLOWSetInitialGuessArray_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, + double *x0_dev) { + PetscErrorCode ierr; + double *x; + auto &resmgr = umpire::ResourceManager::getInstance(); + + PetscFunctionBegin; + // ierr = PetscPrintf(MPI_COMM_SELF,"Entered + // OPFLOWInitialization\n");CHKERRQ(ierr); + // Do initialization on Host + ierr = (*opflow->modelops.setinitialguess)(opflow, opflow->X, opflow->Lambda); + CHKERRQ(ierr); + ierr = VecGetArray(opflow->X, &x); + CHKERRQ(ierr); + + // Copy from host to device + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + registerWith(x, opflow->nx, resmgr, h_allocator_); + resmgr.copy(x0_dev, x); + + ierr = VecRestoreArray(opflow->X, &x); + CHKERRQ(ierr); + + // ierr = PetscPrintf(MPI_COMM_SELF,"Exit Initialization\n");CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +PetscErrorCode +OPFLOWSetVariableBoundsArray_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, double *xl_dev, + double *xu_dev) { + PetscErrorCode ierr; + double *xl, *xu; + auto &resmgr = umpire::ResourceManager::getInstance(); + + PetscFunctionBegin; + + // Compute variable bounds on the host + ierr = (*opflow->modelops.setvariablebounds)(opflow, opflow->Xl, opflow->Xu); + CHKERRQ(ierr); + ierr = VecGetArray(opflow->Xl, &xl); + CHKERRQ(ierr); + ierr = VecGetArray(opflow->Xu, &xu); + CHKERRQ(ierr); + + // Copy host to device + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + registerWith(xl, opflow->nx, resmgr, h_allocator_); + registerWith(xu, opflow->nx, resmgr, h_allocator_); + resmgr.copy(xl_dev, xl); + resmgr.copy(xu_dev, xu); + + ierr = VecRestoreArray(opflow->Xl, &xl); + CHKERRQ(ierr); + ierr = VecRestoreArray(opflow->Xu, &xu); + CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSetConstraintBoundsArray_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, double *gl_dev, double *gu_dev) { + + PetscErrorCode ierr; + double *gl, *gu; + auto &resmgr = umpire::ResourceManager::getInstance(); + + PetscFunctionBegin; + + // Calculate constraint bounds on host + ierr = + (*opflow->modelops.setconstraintbounds)(opflow, opflow->Gl, opflow->Gu); + CHKERRQ(ierr); + + ierr = VecGetArray(opflow->Gl, &gl); + CHKERRQ(ierr); + ierr = VecGetArray(opflow->Gu, &gu); + CHKERRQ(ierr); + + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + registerWith(gl, opflow->ncon, resmgr, h_allocator_); + registerWith(gu, opflow->ncon, resmgr, h_allocator_); + + // Copy from host to device + resmgr.copy(gl_dev, gl); + resmgr.copy(gu_dev, gu); + + ierr = VecRestoreArray(opflow->Gl, &gl); + CHKERRQ(ierr); + ierr = VecRestoreArray(opflow->Gu, &gu); + CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +/** EQUALITY CONSTRAINTS */ +PetscErrorCode OPFLOWComputeEqualityConstraintsArray_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, double *ge_dev) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + BUSParamsRajaHiop *busparams = &pbpolrajahiopsparse->busparams; + GENParamsRajaHiop *genparams = &pbpolrajahiopsparse->genparams; + LOADParamsRajaHiop *loadparams = &pbpolrajahiopsparse->loadparams; + LINEParamsRajaHiop *lineparams = &pbpolrajahiopsparse->lineparams; + PetscErrorCode ierr; + PetscInt flps = 0; + int include_loadloss_variables = opflow->include_loadloss_variables; + int include_powerimbalance_variables = + opflow->include_powerimbalance_variables; + + PetscFunctionBegin; + // PetscPrintf(MPI_COMM_SELF,"Entered Equality constraints\n"); + + // Zero out array + auto &resmgr = umpire::ResourceManager::getInstance(); + resmgr.memset(ge_dev, 0, opflow->nconeq * sizeof(double)); + + /* Generator contributions */ + int *g_gidxbus = genparams->gidxbus_dev_; + int *g_xidx = genparams->xidx_dev_; + RAJA::forall( + RAJA::RangeSegment(0, genparams->ngenON), + RAJA_LAMBDA(RAJA::Index_type i) { + RAJA::atomicSub(&ge_dev[g_gidxbus[i]], + x_dev[g_xidx[i]]); + RAJA::atomicSub(&ge_dev[g_gidxbus[i] + 1], + x_dev[g_xidx[i] + 1]); + }); + flps += genparams->ngenON * 2; + + if (opflow->has_gensetpoint) { + int *g_geqidxgen = genparams->geqidxgen_dev_; + double *g_pgs = genparams->pgs_dev_; + RAJA::forall( + RAJA::RangeSegment(0, genparams->ngenON), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pg, delPg, Pgset; + Pg = x_dev[g_xidx[i]]; + delPg = x_dev[g_xidx[i] + 2]; + Pgset = x_dev[g_xidx[i] + 3]; + + ge_dev[g_geqidxgen[i]] = Pgset + delPg - Pg; + ge_dev[g_geqidxgen[i] + 1] = Pgset - g_pgs[i]; + }); + } + + /* Load contributions */ + double *pl = loadparams->pl_dev_; + double *ql = loadparams->ql_dev_; + int *l_gidx = loadparams->gidx_dev_; + int *l_xidx = loadparams->xidx_dev_; + RAJA::forall( + RAJA::RangeSegment(0, loadparams->nload), + RAJA_LAMBDA(RAJA::Index_type i) { + if (include_loadloss_variables) { + RAJA::atomicAdd(&ge_dev[l_gidx[i]], pl[i]); + RAJA::atomicAdd(&ge_dev[l_gidx[i] + 1], ql[i]); + RAJA::atomicSub(&ge_dev[l_gidx[i]], + x_dev[l_xidx[i]]); + RAJA::atomicSub(&ge_dev[l_gidx[i] + 1], + x_dev[l_xidx[i] + 1]); + } else { + RAJA::atomicAdd(&ge_dev[l_gidx[i]], pl[i]); + RAJA::atomicAdd(&ge_dev[l_gidx[i] + 1], ql[i]); + } + }); + flps += loadparams->nload * 2; + + /* Bus contributions */ + int *isisolated = busparams->isisolated_dev_; + int *ispvpq = busparams->ispvpq_dev_; + double *gl = busparams->gl_dev_; + double *bl = busparams->bl_dev_; + double *vm = busparams->vm_dev_; + double *va = busparams->va_dev_; + int *b_xidx = busparams->xidx_dev_; + int *b_gidx = busparams->gidx_dev_; + int *b_xidxpimb = busparams->xidxpimb_dev_; + + RAJA::forall( + RAJA::RangeSegment(0, busparams->nbus), RAJA_LAMBDA(RAJA::Index_type i) { + double theta = x_dev[b_xidx[i]]; + double Vm = x_dev[b_xidx[i] + 1]; + RAJA::atomicAdd( + &ge_dev[b_gidx[i]], + isisolated[i] * (theta - va[i] * PETSC_PI / 180.0) + + ispvpq[i] * Vm * Vm * gl[i]); + + RAJA::atomicAdd(&ge_dev[b_gidx[i] + 1], + isisolated[i] * (Vm - vm[i]) - + ispvpq[i] * Vm * Vm * bl[i]); + if (include_powerimbalance_variables) { + double Pimb = x_dev[b_xidxpimb[i]]; + double Qimb = x_dev[b_xidxpimb[i] + 1]; + RAJA::atomicAdd(&ge_dev[b_gidx[i]], Pimb); + RAJA::atomicAdd(&ge_dev[b_gidx[i] + 1], Qimb); + } + }); + flps += busparams->nbus * 14; + + /* Line contributions */ + double *Gff = lineparams->Gff_dev_; + double *Gtt = lineparams->Gtt_dev_; + double *Gft = lineparams->Gft_dev_; + double *Gtf = lineparams->Gtf_dev_; + + double *Bff = lineparams->Bff_dev_; + double *Btt = lineparams->Btt_dev_; + double *Bft = lineparams->Bft_dev_; + double *Btf = lineparams->Btf_dev_; + + int *xidxf = lineparams->xidxf_dev_; + int *xidxt = lineparams->xidxt_dev_; + int *geqidxf = lineparams->geqidxf_dev_; + int *geqidxt = lineparams->geqidxt_dev_; + + RAJA::forall( + RAJA::RangeSegment(0, lineparams->nlineON), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pf, Qf, Pt, Qt; + double thetaf = x_dev[xidxf[i]], Vmf = x_dev[xidxf[i] + 1]; + double thetat = x_dev[xidxt[i]], Vmt = x_dev[xidxt[i] + 1]; + double thetaft = thetaf - thetat; + double thetatf = thetat - thetaf; + + Pf = Gff[i] * Vmf * Vmf + + Vmf * Vmt * (Gft[i] * cos(thetaft) + Bft[i] * sin(thetaft)); + Qf = -Bff[i] * Vmf * Vmf + + Vmf * Vmt * (-Bft[i] * cos(thetaft) + Gft[i] * sin(thetaft)); + Pt = Gtt[i] * Vmt * Vmt + + Vmt * Vmf * (Gtf[i] * cos(thetatf) + Btf[i] * sin(thetatf)); + Qt = -Btt[i] * Vmt * Vmt + + Vmt * Vmf * (-Btf[i] * cos(thetatf) + Gtf[i] * sin(thetatf)); + + RAJA::atomicAdd(&ge_dev[geqidxf[i]], Pf); + RAJA::atomicAdd(&ge_dev[geqidxf[i] + 1], Qf); + RAJA::atomicAdd(&ge_dev[geqidxt[i]], Pt); + RAJA::atomicAdd(&ge_dev[geqidxt[i] + 1], Qt); + }); + flps += lineparams->nlineON * + (46 + (4 * EXAGO_FLOPS_COSOP) + (4 * EXAGO_FLOPS_SINOP)); + // PetscPrintf(MPI_COMM_SELF,"Exit Equality Constraints\n"); + + ierr = PetscLogFlops(flps); + CHKERRQ(ierr); + PetscFunctionReturn(0); +} + +/** INEQUALITY CONSTRAINTS **/ +PetscErrorCode OPFLOWComputeInequalityConstraintsArray_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, double *gi_dev) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + LINEParamsRajaHiop *lineparams = &pbpolrajahiopsparse->lineparams; + PetscInt flps = 0; + PetscErrorCode ierr; + + PetscFunctionBegin; + // PetscPrintf(MPI_COMM_SELF,"Entered Inequality Constraints\n"); + + // Zero out array + auto &resmgr = umpire::ResourceManager::getInstance(); + resmgr.memset(gi_dev, 0, opflow->nconineq * sizeof(double)); + + /* Line contributions */ + double *Gff = lineparams->Gff_dev_; + double *Gtt = lineparams->Gtt_dev_; + double *Gft = lineparams->Gft_dev_; + double *Gtf = lineparams->Gtf_dev_; + + double *Bff = lineparams->Bff_dev_; + double *Btt = lineparams->Btt_dev_; + double *Bft = lineparams->Bft_dev_; + double *Btf = lineparams->Btf_dev_; + + int *linelimidx = lineparams->linelimidx_dev_; + int *xidxf = lineparams->xidxf_dev_; + int *xidxt = lineparams->xidxt_dev_; + int *gineqidx = lineparams->gineqidx_dev_; + if (lineparams->nlinelim) { + RAJA::forall( + RAJA::RangeSegment(0, lineparams->nlinelim), + RAJA_LAMBDA(RAJA::Index_type i) { + int j = linelimidx[i]; + double Pf, Qf, Pt, Qt, Sf2, St2; + double thetaf = x_dev[xidxf[j]], Vmf = x_dev[xidxf[j] + 1]; + double thetat = x_dev[xidxt[j]], Vmt = x_dev[xidxt[j] + 1]; + double thetaft = thetaf - thetat; + double thetatf = thetat - thetaf; + + Pf = Gff[j] * Vmf * Vmf + + Vmf * Vmt * (Gft[j] * cos(thetaft) + Bft[j] * sin(thetaft)); + Qf = -Bff[j] * Vmf * Vmf + + Vmf * Vmt * (-Bft[j] * cos(thetaft) + Gft[j] * sin(thetaft)); + Pt = Gtt[j] * Vmt * Vmt + + Vmt * Vmf * (Gtf[j] * cos(thetatf) + Btf[j] * sin(thetatf)); + Qt = -Btt[j] * Vmt * Vmt + + Vmt * Vmf * (-Btf[j] * cos(thetatf) + Gtf[j] * sin(thetatf)); + + Sf2 = Pf * Pf + Qf * Qf; + St2 = Pt * Pt + Qt * Qt; + + gi_dev[gineqidx[i]] = Sf2; + gi_dev[gineqidx[i] + 1] = St2; + }); + flps += lineparams->nlinelim * + (72 + (4 * EXAGO_FLOPS_COSOP) + (4 * EXAGO_FLOPS_SINOP)); + } + // PetscPrintf(MPI_COMM_SELF,"Exit Inequality Constraints\n"); + ierr = PetscLogFlops(flps); + CHKERRQ(ierr); + PetscFunctionReturn(0); +} + +/** OBJECTIVE FUNCTION **/ +// Note: This kernel (and all the kernels for this model assume that the data +// has been already allocated on the device. x_dev is pointer to array on the +// GPU +PetscErrorCode OPFLOWComputeObjectiveArray_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, double *obj) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + GENParamsRajaHiop *genparams = &pbpolrajahiopsparse->genparams; + LOADParamsRajaHiop *loadparams = &pbpolrajahiopsparse->loadparams; + BUSParamsRajaHiop *busparams = &pbpolrajahiopsparse->busparams; + PetscErrorCode ierr; + PS ps = opflow->ps; + int isobj_gencost = opflow->obj_gencost; + double MVAbase = ps->MVAbase; + + PetscFunctionBegin; + // PetscPrintf(MPI_COMM_SELF,"Entered objective function\n"); + + // You need local copies of device pointers so that lambda can capture them + double *cost_alpha = genparams->cost_alpha_dev_; + double *cost_beta = genparams->cost_beta_dev_; + double *cost_gamma = genparams->cost_gamma_dev_; + int *xidx = genparams->xidx_dev_; + int *l_xidx = loadparams->xidx_dev_; + int *b_xidxpimb = busparams->xidxpimb_dev_; + + /* Generator objective function contributions */ + // Set up reduce sum object + RAJA::ReduceSum obj_val_sum(0.0); + // Compute reduction on CUDA device + RAJA::forall( + RAJA::RangeSegment(0, genparams->ngenON), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pg = x_dev[xidx[i]] * MVAbase; + obj_val_sum += isobj_gencost * (cost_alpha[i] * Pg * Pg + + cost_beta[i] * Pg + cost_gamma[i]); + }); + + if (opflow->include_loadloss_variables) { + RAJA::forall( + RAJA::RangeSegment(0, loadparams->nload), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pdloss = x_dev[l_xidx[i]]; + double Qdloss = x_dev[l_xidx[i] + 1]; + obj_val_sum += loadparams->loadloss_penalty_dev_[i] * ps->MVAbase * + ps->MVAbase * (Pdloss * Pdloss + Qdloss * Qdloss); + }); + } + + /* Powerimbalance contributions */ + if (opflow->include_powerimbalance_variables) { + RAJA::forall( + RAJA::RangeSegment(0, busparams->nbus), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pimb = x_dev[b_xidxpimb[i]]; + double Qimb = x_dev[b_xidxpimb[i] + 1]; + obj_val_sum += busparams->powerimbalance_penalty_dev_[i] * + ps->MVAbase * ps->MVAbase * + (Pimb * Pimb + Qimb * Qimb); + }); + } + + *obj = static_cast(obj_val_sum.get()); + ierr = PetscLogFlops(genparams->ngenON * 8.0); + CHKERRQ(ierr); + + // PetscPrintf(MPI_COMM_SELF,"Exit objective function\n"); + PetscFunctionReturn(0); +} + +/** GRADIENT **/ +PetscErrorCode OPFLOWComputeGradientArray_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, double *grad_dev) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + GENParamsRajaHiop *genparams = &pbpolrajahiopsparse->genparams; + LOADParamsRajaHiop *loadparams = &pbpolrajahiopsparse->loadparams; + BUSParamsRajaHiop *busparams = &pbpolrajahiopsparse->busparams; + PS ps = opflow->ps; + int isobj_gencost = opflow->obj_gencost; + double MVAbase = ps->MVAbase; + PetscErrorCode ierr; + + PetscFunctionBegin; + // PetscPrintf(MPI_COMM_SELF,"Entered gradient function\n"); + // Zero out array + auto &resmgr = umpire::ResourceManager::getInstance(); + resmgr.memset(grad_dev, 0); + + /* Generator gradient contributions */ + double *cost_alpha = genparams->cost_alpha_dev_; + double *cost_beta = genparams->cost_beta_dev_; + int *xidx = genparams->xidx_dev_; + int *l_xidx = loadparams->xidx_dev_; + int *b_xidxpimb = busparams->xidxpimb_dev_; + + RAJA::forall( + RAJA::RangeSegment(0, genparams->ngenON), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pg = x_dev[xidx[i]] * MVAbase; + grad_dev[xidx[i]] = + isobj_gencost * MVAbase * (2.0 * cost_alpha[i] * Pg + cost_beta[i]); + }); + + if (opflow->include_loadloss_variables) { + RAJA::forall( + RAJA::RangeSegment(0, loadparams->nload), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pdloss = x_dev[l_xidx[i]]; + double Qdloss = x_dev[l_xidx[i] + 1]; + grad_dev[l_xidx[i]] = loadparams->loadloss_penalty_dev_[i] * + ps->MVAbase * ps->MVAbase * 2 * Pdloss; + grad_dev[l_xidx[i] + 1] = loadparams->loadloss_penalty_dev_[i] * + ps->MVAbase * ps->MVAbase * 2 * Qdloss; + }); + } + + if (opflow->include_powerimbalance_variables) { + RAJA::forall( + RAJA::RangeSegment(0, busparams->nbus), + RAJA_LAMBDA(RAJA::Index_type i) { + double Pimb = x_dev[b_xidxpimb[i]]; + double Qimb = x_dev[b_xidxpimb[i] + 1]; + grad_dev[b_xidxpimb[i]] = busparams->powerimbalance_penalty_dev_[i] * + ps->MVAbase * ps->MVAbase * 2 * Pimb; + grad_dev[b_xidxpimb[i] + 1] = + busparams->powerimbalance_penalty_dev_[i] * ps->MVAbase * + ps->MVAbase * 2 * Qimb; + }); + } + + ierr = PetscLogFlops(genparams->ngenON * 6.0); + CHKERRQ(ierr); + // PetscPrintf(MPI_COMM_SELF,"Exit gradient function\n"); + + PetscFunctionReturn(0); +} + +PetscErrorCode +OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, int *iJacS_dev, int *jJacS_dev, + double *MJacS_dev) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + PetscErrorCode ierr; + double *x, *values; + PetscInt *iRowstart, *jColstart; + PetscInt roffset, coffset; + PetscInt nrow, ncol; + PetscInt nvals; + const PetscInt *cols; + const PetscScalar *vals; + PetscInt i, j; + auto &resmgr = umpire::ResourceManager::getInstance(); + + PetscFunctionBegin; + + if (MJacS_dev == NULL) { + /* Set locations only */ + + if (opflow->Nconineq) { + // Create arrays on host to store i,j, and val arrays + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + + pbpolrajahiopsparse->i_jacineq = + (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); + pbpolrajahiopsparse->j_jacineq = + (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); + pbpolrajahiopsparse->val_jacineq = (double *)(h_allocator_.allocate( + opflow->nnz_ineqjacsp * sizeof(double))); + + iRowstart = pbpolrajahiopsparse->i_jacineq; + jColstart = pbpolrajahiopsparse->j_jacineq; + + /* Inequality constraints start after equality constraints + Hence the offset + */ + roffset = opflow->nconeq; + coffset = 0; + + ierr = (*opflow->modelops.computeinequalityconstraintjacobian)( + opflow, opflow->X, opflow->Jac_Gi); + CHKERRQ(ierr); + + ierr = MatGetSize(opflow->Jac_Gi, &nrow, &ncol); + CHKERRQ(ierr); + /* Copy over locations to triplet format */ + for (i = 0; i < nrow; i++) { + ierr = MatGetRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + for (j = 0; j < nvals; j++) { + iRowstart[j] = roffset + i; + jColstart[j] = coffset + cols[j]; + } + /* Increment iRow,jCol pointers */ + iRowstart += nvals; + jColstart += nvals; + ierr = MatRestoreRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + } + + // Copy over i_jacineq and j_jacineq arrays to device + resmgr.copy(iJacS_dev + opflow->nnz_eqjacsp, + pbpolrajahiopsparse->i_jacineq); + resmgr.copy(jJacS_dev + opflow->nnz_eqjacsp, + pbpolrajahiopsparse->j_jacineq); + + ierr = PetscLogEventEnd(opflow->ineqconsjaclogger, 0, 0, 0, 0); + CHKERRQ(ierr); + } + } else { + if (opflow->Nconineq) { + ierr = PetscLogEventBegin(opflow->ineqconsjaclogger, 0, 0, 0, 0); + CHKERRQ(ierr); + + ierr = VecGetArray(opflow->X, &x); + CHKERRQ(ierr); + + // Copy from device to host + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + registerWith(x, opflow->nx, resmgr, h_allocator_); + resmgr.copy((double *)x, (double *)x_dev); + + ierr = VecRestoreArray(opflow->X, &x); + CHKERRQ(ierr); + + /* Compute inequality constraint jacobian */ + ierr = (*opflow->modelops.computeinequalityconstraintjacobian)( + opflow, opflow->X, opflow->Jac_Gi); + CHKERRQ(ierr); + + ierr = MatGetSize(opflow->Jac_Gi, &nrow, &ncol); + CHKERRQ(ierr); + + values = pbpolrajahiopsparse->val_jacineq; + /* Copy over values */ + for (i = 0; i < nrow; i++) { + ierr = MatGetRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + for (j = 0; j < nvals; j++) { + values[j] = vals[j]; + } + values += nvals; + ierr = MatRestoreRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + } + // Copy over val_jacineq to device + resmgr.copy(MJacS_dev + opflow->nnz_eqjacsp, + pbpolrajahiopsparse->val_jacineq); + + ierr = PetscLogEventEnd(opflow->ineqconsjaclogger, 0, 0, 0, 0); + CHKERRQ(ierr); + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode +OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, int *iJacS_dev, int *jJacS_dev, + double *MJacS_dev) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + PetscErrorCode ierr; + PetscInt *iRowstart, *jColstart; + PetscScalar *x, *values; + PetscInt roffset, coffset; + PetscInt nrow, ncol; + PetscInt nvals; + const PetscInt *cols; + const PetscScalar *vals; + PetscInt i, j; + auto &resmgr = umpire::ResourceManager::getInstance(); + + PetscFunctionBegin; + + if (MJacS_dev == NULL) { + /* Set locations only */ + + roffset = 0; + coffset = 0; + + // Create arrays on host to store i,j, and val arrays + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + + pbpolrajahiopsparse->i_jaceq = + (int *)(h_allocator_.allocate(opflow->nnz_eqjacsp * sizeof(int))); + pbpolrajahiopsparse->j_jaceq = + (int *)(h_allocator_.allocate(opflow->nnz_eqjacsp * sizeof(int))); + pbpolrajahiopsparse->val_jaceq = + (double *)(h_allocator_.allocate(opflow->nnz_eqjacsp * sizeof(double))); + + iRowstart = pbpolrajahiopsparse->i_jaceq; + jColstart = pbpolrajahiopsparse->j_jaceq; + + ierr = (*opflow->modelops.computeequalityconstraintjacobian)( + opflow, opflow->X, opflow->Jac_Ge); + CHKERRQ(ierr); + + ierr = MatGetSize(opflow->Jac_Ge, &nrow, &ncol); + CHKERRQ(ierr); + + /* Copy over locations to triplet format */ + for (i = 0; i < nrow; i++) { + ierr = MatGetRow(opflow->Jac_Ge, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + for (j = 0; j < nvals; j++) { + iRowstart[j] = roffset + i; + jColstart[j] = coffset + cols[j]; + } + /* Increment iRow,jCol pointers */ + iRowstart += nvals; + jColstart += nvals; + ierr = MatRestoreRow(opflow->Jac_Ge, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + } + + // Copy over i_jaceq and j_jaceq arrays to device + resmgr.copy(iJacS_dev, pbpolrajahiopsparse->i_jaceq); + resmgr.copy(jJacS_dev, pbpolrajahiopsparse->j_jaceq); + } else { + ierr = PetscLogEventBegin(opflow->eqconsjaclogger, 0, 0, 0, 0); + CHKERRQ(ierr); + + ierr = VecGetArray(opflow->X, &x); + CHKERRQ(ierr); + + // Copy from device to host + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + registerWith(x, opflow->nx, resmgr, h_allocator_); + resmgr.copy((double *)x, (double *)x_dev); + + ierr = VecRestoreArray(opflow->X, &x); + CHKERRQ(ierr); + + /* Compute equality constraint jacobian */ + ierr = (*opflow->modelops.computeequalityconstraintjacobian)( + opflow, opflow->X, opflow->Jac_Ge); + CHKERRQ(ierr); + + ierr = MatGetSize(opflow->Jac_Ge, &nrow, &ncol); + CHKERRQ(ierr); + + values = pbpolrajahiopsparse->val_jaceq; + + /* Copy over values */ + for (i = 0; i < nrow; i++) { + ierr = MatGetRow(opflow->Jac_Ge, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + for (j = 0; j < nvals; j++) { + values[j] = vals[j]; + } + values += nvals; + ierr = MatRestoreRow(opflow->Jac_Ge, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + } + + // Copy over val_ineq to device + resmgr.copy(MJacS_dev, pbpolrajahiopsparse->val_jaceq); + + ierr = PetscLogEventEnd(opflow->eqconsjaclogger, 0, 0, 0, 0); + CHKERRQ(ierr); + } + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWComputeSparseHessian_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *x_dev, const double *lambda_dev, int *iHSS_dev, + int *jHSS_dev, double *MHSS_dev) { + PbpolModelRajaHiop *pbpolrajahiopsparse = + reinterpret_cast(opflow->model); + PetscErrorCode ierr; + PetscInt *iRow, *jCol; + PetscScalar *x, *values, *lambda; + PetscInt nrow; + PetscInt nvals; + const PetscInt *cols; + const PetscScalar *vals; + PetscInt i, j; + PetscInt ctr = 0; + auto &resmgr = umpire::ResourceManager::getInstance(); + + PetscFunctionBegin; + + if (iHSS_dev != NULL && jHSS_dev != NULL) { + + // Create arrays on host to store i,j, and val arrays + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + + pbpolrajahiopsparse->i_hess = + (int *)(h_allocator_.allocate(opflow->nnz_hesssp * sizeof(int))); + pbpolrajahiopsparse->j_hess = + (int *)(h_allocator_.allocate(opflow->nnz_hesssp * sizeof(int))); + pbpolrajahiopsparse->val_hess = + (double *)(h_allocator_.allocate(opflow->nnz_hesssp * sizeof(double))); + + iRow = pbpolrajahiopsparse->i_hess; + jCol = pbpolrajahiopsparse->j_hess; + + ierr = (*opflow->modelops.computehessian)( + opflow, opflow->X, opflow->Lambdae, opflow->Lambdai, opflow->Hes); + CHKERRQ(ierr); + ierr = MatGetSize(opflow->Hes, &nrow, &nrow); + CHKERRQ(ierr); + + /* Copy over locations to triplet format */ + /* Note that HIOP requires a upper triangular Hessian as oppposed + to IPOPT which requires a lower triangular Hessian + */ + for (i = 0; i < nrow; i++) { + ierr = MatGetRow(opflow->Hes, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + ctr = 0; + for (j = 0; j < nvals; j++) { + if (cols[j] >= i) { /* upper triangle */ + /* save as upper triangle locations */ + iRow[ctr] = i; + jCol[ctr] = cols[j]; + ctr++; + } + } + iRow += ctr; + jCol += ctr; + ierr = MatRestoreRow(opflow->Hes, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + } + + // Copy over i_hess and j_hess arrays to device + resmgr.copy(iHSS_dev, pbpolrajahiopsparse->i_hess); + resmgr.copy(jHSS_dev, pbpolrajahiopsparse->j_hess); + } else { + + ierr = VecGetArray(opflow->X, &x); + CHKERRQ(ierr); + + // Copy from device to host + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + registerWith(x, opflow->nx, resmgr, h_allocator_); + resmgr.copy((double *)x, (double *)x_dev); + + ierr = VecRestoreArray(opflow->X, &x); + CHKERRQ(ierr); + + ierr = VecGetArray(opflow->Lambda, &lambda); + + registerWith(lambda, opflow->ncon, resmgr, h_allocator_); + // copy lambda from device to host + resmgr.copy((double *)lambda, (double *)lambda_dev); + + ierr = VecPlaceArray(opflow->Lambdae, lambda); + CHKERRQ(ierr); + if (opflow->Nconineq) { + ierr = VecPlaceArray(opflow->Lambdai, lambda + opflow->nconeq); + CHKERRQ(ierr); + } + + /* Compute Hessian */ + ierr = (*opflow->modelops.computehessian)( + opflow, opflow->X, opflow->Lambdae, opflow->Lambdai, opflow->Hes); + CHKERRQ(ierr); + + ierr = VecResetArray(opflow->Lambdae); + CHKERRQ(ierr); + if (opflow->Nconineq) { + ierr = VecResetArray(opflow->Lambdai); + CHKERRQ(ierr); + } + + ierr = VecRestoreArray(opflow->Lambda, &lambda); + CHKERRQ(ierr); + + if (opflow->modelops.computeauxhessian) { + ierr = VecGetArray(opflow->X, &x); + CHKERRQ(ierr); + ierr = (*opflow->modelops.computeauxhessian)(opflow, x, opflow->Hes, + opflow->userctx); + CHKERRQ(ierr); + ierr = VecRestoreArray(opflow->X, &x); + CHKERRQ(ierr); + + ierr = MatAssemblyBegin(opflow->Hes, MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); + ierr = MatAssemblyEnd(opflow->Hes, MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); + } + + /* Copy over values */ + ierr = MatGetSize(opflow->Hes, &nrow, &nrow); + CHKERRQ(ierr); + + values = pbpolrajahiopsparse->val_hess; + + for (i = 0; i < nrow; i++) { + ierr = MatGetRow(opflow->Hes, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + ctr = 0; + for (j = 0; j < nvals; j++) { + if (cols[j] >= i) { /* Upper triangle values (same as lower triangle) */ + values[ctr] = vals[j]; + ctr++; + } + } + values += ctr; + ierr = MatRestoreRow(opflow->Hes, i, &nvals, &cols, &vals); + CHKERRQ(ierr); + } + + // Copy over val_ineq to device + resmgr.copy(MHSS_dev, pbpolrajahiopsparse->val_hess); + } + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolutionCallback_PBPOLRAJAHIOPSPARSE( + OPFLOW opflow, const double *xsol, const double *z_L, const double *z_U, + const double *gsol, const double *lamsol, double obj_value) { + (void)z_L; + (void)z_U; + (void)obj_value; + + PetscErrorCode ierr; + PetscScalar *x, *lam, *g; + + auto &resmgr = umpire::ResourceManager::getInstance(); + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + + /* FIXME: It is assumed that array arguments are on the device. If + there are other posibilities, opflow->mem_space should be set + properly and used to decide how to get at those arrays */ + + ierr = VecGetArray(opflow->X, &x); + CHKERRQ(ierr); + /* Copy xsol from device to host */ + resmgr.copy(x, (double *)xsol); + ierr = VecRestoreArray(opflow->X, &x); + CHKERRQ(ierr); + + if (lamsol) { + /* HIOP returns a NULL for lamsol - probably lamsol needs to be added to + HIOP. Need to remove this condition once it is fixed + */ + ierr = VecGetArray(opflow->Lambda, &lam); + CHKERRQ(ierr); + + /* Create temporary vectors for copying values to HOST */ + double *lamsol_host = + (double *)h_allocator_.allocate(opflow->ncon * sizeof(double)); + + /* Copy lamsol from device to host */ + resmgr.copy(lamsol_host, (double *)lamsol); + + ierr = PetscMemcpy(lam, (double *)lamsol_host, + opflow->nconeq * sizeof(PetscScalar)); + CHKERRQ(ierr); + if (opflow->Nconineq) { + ierr = PetscMemcpy(lam + opflow->nconeq, + (double *)(lamsol_host + opflow->nconeq), + opflow->nconineq * sizeof(PetscScalar)); + CHKERRQ(ierr); + } + ierr = VecRestoreArray(opflow->Lambda, &lam); + CHKERRQ(ierr); + h_allocator_.deallocate(lamsol_host); + } else { + ierr = VecSet(opflow->Lambda, -9999.0); + CHKERRQ(ierr); + } + + if (gsol) { + ierr = VecGetArray(opflow->G, &g); + CHKERRQ(ierr); + + /* Create temporary vectors for copying values to HOST */ + double *gsol_host = + (double *)h_allocator_.allocate(opflow->ncon * sizeof(double)); + + /* Copy gsol from device to host */ + resmgr.copy(gsol_host, (double *)gsol); + + ierr = PetscMemcpy(g, (double *)gsol_host, + opflow->nconeq * sizeof(PetscScalar)); + CHKERRQ(ierr); + if (opflow->Nconineq) { + ierr = PetscMemcpy(g + opflow->nconeq, + (double *)(gsol_host + opflow->nconeq), + opflow->nconineq * sizeof(PetscScalar)); + CHKERRQ(ierr); + } + ierr = VecRestoreArray(opflow->G, &g); + CHKERRQ(ierr); + h_allocator_.deallocate(gsol_host); + } + PetscFunctionReturn(0); +} + +#endif +#endif diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.hpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.hpp new file mode 100644 index 00000000..9b941f3e --- /dev/null +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.hpp @@ -0,0 +1,32 @@ +#include + +#if defined(EXAGO_ENABLE_RAJA) +#if defined(EXAGO_ENABLE_HIOP_SPARSE) + +#pragma once + +#include +#include + +#include "pbpolrajahiopsparse.hpp" + +/** + + @brief Takes a raw pointer to some data and registers it with the + allocator and resource manager. + + TODO: Ensure that we do not have to deregister the allocation with + the resource manager. This is done with `resmgr.deregisterAllocation(T* ptr)`, + but must we do this with ever allocation? Will we incurr a memory leak if not? + +*/ +template +void registerWith(T *ptr, SizeType N, umpire::ResourceManager &resmgr, + umpire::Allocator &allocator) { + umpire::util::AllocationRecord record{ptr, sizeof(T) * N, + allocator.getAllocationStrategy()}; + resmgr.registerAllocation(ptr, record); +} + +#endif // EXAGO_ENABLE_HIOP_SPARSE +#endif // EXAGO_ENABLE_HIOP_RAJA diff --git a/src/opflow/solver/hiop/opflow_hiopsparse.cpp b/src/opflow/solver/hiop/opflow_hiopsparse.cpp index eed864f3..5b9ca60b 100644 --- a/src/opflow/solver/hiop/opflow_hiopsparse.cpp +++ b/src/opflow/solver/hiop/opflow_hiopsparse.cpp @@ -753,6 +753,7 @@ PetscErrorCode OPFLOWSolverDestroy_HIOPSPARSE(OPFLOW opflow) { PetscFunctionBegin; + delete hiop->solver; delete hiop->sp; delete hiop->nlp; diff --git a/src/opflow/solver/hiop/opflow_hiopsparse.h b/src/opflow/solver/hiop/opflow_hiopsparse.h index e86b09e9..80e3e2c1 100644 --- a/src/opflow/solver/hiop/opflow_hiopsparse.h +++ b/src/opflow/solver/hiop/opflow_hiopsparse.h @@ -16,10 +16,6 @@ #include #endif -#ifdef HIOP_USE_MAGMA -#include "magma_v2.h" -#endif - class OPFLOWHIOPSPARSEInterface : public hiop::hiopInterfaceSparse { public: OPFLOWHIOPSPARSEInterface(OPFLOW); diff --git a/src/opflow/solver/hiop/opflow_hiopsparsegpu.cpp b/src/opflow/solver/hiop/opflow_hiopsparsegpu.cpp new file mode 100644 index 00000000..3517c98a --- /dev/null +++ b/src/opflow/solver/hiop/opflow_hiopsparsegpu.cpp @@ -0,0 +1,527 @@ +#include +#if defined(EXAGO_ENABLE_HIOP) +#if defined(EXAGO_ENABLE_HIOP_SPARSE) + +#include +#include "opflow_hiopsparsegpu.hpp" + +OPFLOWHIOPSPARSEGPUInterface::OPFLOWHIOPSPARSEGPUInterface(OPFLOW opflowin) { + opflow = opflowin; +} + +bool OPFLOWHIOPSPARSEGPUInterface::get_prob_sizes(hiop::size_type &n, + hiop::size_type &m) { + n = opflow->nx; + m = opflow->ncon; + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::get_vars_info(const hiop::size_type &n, + double *xlow, double *xupp, + NonlinearityType *type) { + PetscErrorCode ierr; + PetscInt i; + + ierr = (*opflow->modelops.setvariableboundsarray)(opflow, xlow, xupp); + CHKERRQ(ierr); + + for (i = 0; i < n; i++) { + type[i] = hiopNonlinear; + } + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::get_cons_info(const hiop::size_type &m, + double *clow, double *cupp, + NonlinearityType *type) { + PetscInt i; + PetscErrorCode ierr; + + ierr = (*opflow->modelops.setconstraintboundsarray)(opflow, clow, cupp); + CHKERRQ(ierr); + + for (i = 0; i < m; i++) + type[i] = hiopNonlinear; + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::get_sparse_blocks_info( + hiop::size_type &nx, hiop::size_type &nnz_sparse_Jaceq, + hiop::size_type &nnz_sparse_Jacineq, + hiop::size_type &nnz_sparse_Hess_Lagr) { + PetscErrorCode ierr; + MatInfo info_eq, info_ineq, info_hes; + + nx = opflow->nx; + + /* Compute nonzeros for the Jacobian */ + /* Equality constraint Jacobian */ + ierr = (*opflow->modelops.computeequalityconstraintjacobian)( + opflow, opflow->X, opflow->Jac_Ge); + CHKERRQ(ierr); + ierr = MatSetOption(opflow->Jac_Ge, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE); + CHKERRQ(ierr); + + ierr = MatGetInfo(opflow->Jac_Ge, MAT_LOCAL, &info_eq); + CHKERRQ(ierr); + + nnz_sparse_Jaceq = opflow->nnz_eqjacsp = info_eq.nz_used; + + nnz_sparse_Jacineq = 0; + if (opflow->Nconineq) { + ierr = (*opflow->modelops.computeinequalityconstraintjacobian)( + opflow, opflow->X, opflow->Jac_Gi); + CHKERRQ(ierr); + ierr = + MatSetOption(opflow->Jac_Gi, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE); + CHKERRQ(ierr); + + ierr = MatGetInfo(opflow->Jac_Gi, MAT_LOCAL, &info_ineq); + CHKERRQ(ierr); + + nnz_sparse_Jacineq = opflow->nnz_ineqjacsp = info_ineq.nz_used; + } + + /* Compute non-zeros for Hessian */ + ierr = (*opflow->modelops.computehessian)(opflow, opflow->X, opflow->Lambdae, + opflow->Lambdai, opflow->Hes); + CHKERRQ(ierr); + ierr = MatSetOption(opflow->Hes, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE); + CHKERRQ(ierr); + + ierr = MatGetInfo(opflow->Hes, MAT_LOCAL, &info_hes); + CHKERRQ(ierr); + + nnz_sparse_Hess_Lagr = (info_hes.nz_used - opflow->nx) / 2 + opflow->nx; + + opflow->nnz_hesssp = nnz_sparse_Hess_Lagr; + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_f(const hiop::size_type &n, + const double *x, bool new_x, + double &obj_value) { + PetscErrorCode ierr; + + (void)n; + (void)new_x; + + obj_value = 0.0; + + /* Compute objective */ + ierr = (*opflow->modelops.computeobjectivearray)(opflow, x, &obj_value); + CHKERRQ(ierr); + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_cons(const hiop::size_type &n, + const hiop::size_type &m, + const hiop::size_type &num_cons, + const hiop::size_type *idx_cons, + const double *x, bool new_x, + double *cons) { + (void)n; + (void)m; + (void)num_cons; + (void)idx_cons; + (void)x; + (void)new_x; + (void)cons; + return false; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_cons(const hiop::size_type &n, + const hiop::size_type &m, + const double *x, bool new_x, + double *cons) { + PetscErrorCode ierr; + + (void)n; + (void)m; + (void)new_x; + (void)cons; + + /* Equality constaints */ + ierr = PetscLogEventBegin(opflow->eqconslogger, 0, 0, 0, 0); + CHKERRQ(ierr); + + ierr = (*opflow->modelops.computeequalityconstraintsarray)(opflow, x, cons); + CHKERRQ(ierr); + ierr = PetscLogEventEnd(opflow->eqconslogger, 0, 0, 0, 0); + CHKERRQ(ierr); + + if (opflow->nconineq) { + ierr = PetscLogEventBegin(opflow->ineqconslogger, 0, 0, 0, 0); + CHKERRQ(ierr); + /* Inequality constraints */ + ierr = (*opflow->modelops.computeinequalityconstraintsarray)( + opflow, x, cons + opflow->nconeq); + CHKERRQ(ierr); + ierr = PetscLogEventEnd(opflow->ineqconslogger, 0, 0, 0, 0); + CHKERRQ(ierr); + } + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_grad_f(const hiop::size_type &n, + const double *x, bool new_x, + double *gradf) { + PetscErrorCode ierr; + + (void)n; + (void)new_x; + + ierr = PetscLogEventBegin(opflow->gradlogger, 0, 0, 0, 0); + CHKERRQ(ierr); + ierr = (*opflow->modelops.computegradientarray)(opflow, x, gradf); + CHKERRQ(ierr); + ierr = PetscLogEventEnd(opflow->gradlogger, 0, 0, 0, 0); + CHKERRQ(ierr); + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_Jac_cons( + const hiop::size_type &n, const hiop::size_type &m, + const hiop::size_type &num_cons, const hiop::size_type *idx_cons, + const double *x, bool new_x, const hiop::size_type &nnzJacS, + hiop::index_type *iJacS_dev, hiop::index_type *jJacS_dev, + double *MJacS_dev) { + + (void)n; + (void)m; + (void)num_cons; + (void)idx_cons; + (void)x; + (void)new_x; + (void)nnzJacS; + (void)iJacS_dev; + (void)jJacS_dev; + (void)MJacS_dev; + + return false; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_Jac_cons( + const hiop::size_type &n, const hiop::size_type &m, const double *x, + bool new_x, const hiop::size_type &nnzJacS, hiop::index_type *iJacS_dev, + hiop::index_type *jJacS_dev, double *MJacS_dev) { + PetscErrorCode ierr; + + (void)n; + (void)m; + (void)new_x; + (void)nnzJacS; + + /* Sparse Jacobian */ + ierr = (*opflow->modelops.computesparseequalityconstraintjacobianhiop)( + opflow, x, iJacS_dev, jJacS_dev, MJacS_dev); + CHKERRQ(ierr); + + if (opflow->nconineq) { + /* Sparse Inequality constraint Jacobian */ + ierr = (*opflow->modelops.computesparseinequalityconstraintjacobianhiop)( + opflow, x, iJacS_dev, jJacS_dev, MJacS_dev); + CHKERRQ(ierr); + } + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::eval_Hess_Lagr( + const hiop::size_type &n, const hiop::size_type &m, const double *x_dev, + bool new_x, const double &obj_factor, const double *lambda_dev, + bool new_lambda, const hiop::size_type &nnzHSS, hiop::index_type *iHSS_dev, + hiop::index_type *jHSS_dev, double *MHSS_dev) { + PetscErrorCode ierr; + + (void)n; + (void)m; + (void)new_x; + (void)new_lambda; + (void)nnzHSS; + + opflow->obj_factor = obj_factor; + + /* Compute sparse hessian */ + ierr = PetscLogEventBegin(opflow->sparsehesslogger, 0, 0, 0, 0); + CHKERRQ(ierr); + ierr = (*opflow->modelops.computesparsehessianhiop)( + opflow, x_dev, lambda_dev, iHSS_dev, jHSS_dev, MHSS_dev); + CHKERRQ(ierr); + ierr = PetscLogEventEnd(opflow->sparsehesslogger, 0, 0, 0, 0); + CHKERRQ(ierr); + + return true; +} + +bool OPFLOWHIOPSPARSEGPUInterface::get_starting_point( + const hiop::size_type &global_n, double *x0) { + PetscErrorCode ierr; + + (void)global_n; + + /* Set initial guess */ + ierr = (*opflow->modelops.setinitialguessarray)(opflow, x0); + CHKERRQ(ierr); + + return true; +} + +void OPFLOWHIOPSPARSEGPUInterface::solution_callback( + hiop::hiopSolveStatus status, int n, const double *xsol, const double *z_L, + const double *z_U, int m, const double *gsol, const double *lamsol, + double obj_value) { + (void)n; + (void)m; + PetscErrorCode ierr; + OPFLOWSolver_HIOPSPARSEGPU hiop = (OPFLOWSolver_HIOPSPARSEGPU)opflow->solver; + + /* Copy over solution details */ + hiop->status = status; + opflow->obj = obj_value; + + if (opflow->modelops.solutioncallbackhiop) { + ierr = (*opflow->modelops.solutioncallbackhiop)(opflow, xsol, z_L, z_U, + gsol, lamsol, obj_value); + CHKERRV(ierr); + } +} + +bool OPFLOWHIOPSPARSEGPUInterface::iterate_callback( + int iter, double obj_value, double logbar_obj_value, int n, const double *x, + const double *z_L, const double *z_U, int m_ineq, const double *s, int m, + const double *g, const double *lambda, double inf_pr, double inf_du, + double onenorm_pr_, double mu, double alpha_du, double alpha_pr, + int ls_trials) { + (void)obj_value; + (void)logbar_obj_value; + (void)n; + (void)x; + (void)z_L; + (void)z_U; + (void)m_ineq; + (void)s; + (void)m; + (void)g; + (void)lambda; + (void)inf_pr; + (void)inf_du; + (void)onenorm_pr_; + (void)mu; + (void)alpha_du; + (void)alpha_pr; + (void)ls_trials; + opflow->numits = iter; + return true; +} + +PetscErrorCode OPFLOWSolverSetUp_HIOPSPARSEGPU(OPFLOW opflow) { + PetscErrorCode ierr; + OPFLOWSolver_HIOPSPARSEGPU hiop = (OPFLOWSolver_HIOPSPARSEGPU)opflow->solver; + PetscBool flg1; + int verbose_level = 3; + + PetscFunctionBegin; + + hiop->nlp = new OPFLOWHIOPSPARSEGPUInterface(opflow); + hiop->sp = new hiop::hiopNlpSparse(*hiop->nlp); + + hiop->ipopt_debug = PETSC_FALSE; + + PetscOptionsBegin(opflow->comm->type, NULL, "HIOP options", NULL); + + ierr = PetscOptionsInt("-hiop_verbosity_level", + "HIOP verbosity level (Integer 0 to 12)", "", + verbose_level, &verbose_level, NULL); + CHKERRQ(ierr); +#if defined(EXAGO_ENABLE_IPOPT) + ierr = PetscOptionsBool("-hiop_ipopt_debug", + "Flag enabling debugging HIOP code with IPOPT", "", + hiop->ipopt_debug, &hiop->ipopt_debug, NULL); + CHKERRQ(ierr); +#endif + PetscOptionsEnd(); + +#if defined(EXAGO_ENABLE_IPOPT) + // IPOPT Adapter + if (hiop->ipopt_debug) { + std::cout << "using IPOPT adapter...\n\n"; + hiop->ipoptTNLP = new hiop::hiopSparse2IpoptTNLP(hiop->nlp); + hiop->ipoptApp = new Ipopt::IpoptApplication(); + + // Using options included in HiOp's IpoptAdapter_driver.cpp + hiop->ipoptApp->Options()->SetStringValue("recalc_y", "no"); + hiop->ipoptApp->Options()->SetStringValue("mu_strategy", "monotone"); + hiop->ipoptApp->Options()->SetNumericValue("bound_frac", 1e-8); + hiop->ipoptApp->Options()->SetNumericValue("bound_push", 1e-8); + hiop->ipoptApp->Options()->SetNumericValue("bound_relax_factor", 0.); + hiop->ipoptApp->Options()->SetNumericValue("constr_mult_init_max", 0.001); + hiop->ipoptApp->Options()->SetStringValue("derivative_test", + "second-order"); + + Ipopt::ApplicationReturnStatus status = hiop->ipoptApp->Initialize(); + + if (status != Solve_Succeeded) { + std::cout << std::endl + << std::endl + << "*** Error during initialization!" << std::endl; + return (int)status; + } + PetscFunctionReturn(0); + } +#endif + + hiop->sp->options->SetStringValue("mem_space", "device"); + hiop->sp->options->SetStringValue("compute_mode", "gpu"); + + hiop->sp->options->SetStringValue("dualsInitialization", "zero"); + hiop->sp->options->SetStringValue("duals_init", "zero"); + + hiop->sp->options->SetStringValue("fact_acceptor", "inertia_free"); + hiop->sp->options->SetStringValue("linsol_mode", "speculative"); + hiop->sp->options->SetStringValue("fixed_var", "relax"); + hiop->sp->options->SetStringValue("Hessian", "analytical_exact"); + hiop->sp->options->SetStringValue("KKTLinsys", "xdycyd"); + + hiop->sp->options->SetIntegerValue("verbosity_level", verbose_level); + hiop->sp->options->SetNumericValue("mu0", 1e-1); + hiop->sp->options->SetNumericValue("tolerance", opflow->tolerance); + hiop->sp->options->SetNumericValue("bound_relax_perturb", 1e-4); + hiop->sp->options->SetStringValue("scaling_type", "none"); + + hiop->solver = new hiop::hiopAlgFilterIPMNewton(hiop->sp); + + /* Error if model is not power balance hiop */ + ierr = PetscStrcmp(opflow->modelname.c_str(), OPFLOWMODEL_PBPOLRAJAHIOPSPARSE, + &flg1); + CHKERRQ(ierr); + if (!flg1) { + SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, + "Only PBPOLRAJAHIOPSPARSE model allowed with solver HIOPSPARSEGPU " + "\n Run with -opflow_model " + "PBPOLRAJAHIOPSPARSE\n"); + exit(1); + } + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverSolve_HIOPSPARSEGPU(OPFLOW opflow) { + OPFLOWSolver_HIOPSPARSEGPU hiop = (OPFLOWSolver_HIOPSPARSEGPU)opflow->solver; + + PetscFunctionBegin; +#if defined(EXAGO_ENABLE_IPOPT) + if (!hiop->ipopt_debug) { + hiop->status = hiop->solver->run(); + } else { // Ipopt Adapter + std::cout << "Solving with IPOPT adapter...\n\n"; + ApplicationReturnStatus status = + hiop->ipoptApp->OptimizeTNLP(hiop->ipoptTNLP); + + if (status == Solve_Succeeded) { + std::cout << std::endl + << std::endl + << "*** The problem solved!" << std::endl; + } else { + std::cout << std::endl + << std::endl + << "*** The problem FAILED!" << std::endl; + PetscFunctionReturn(1); + } + } +#else + hiop->status = hiop->solver->run(); +#endif + + PetscFunctionReturn(0); +} + +PetscErrorCode +OPFLOWSolverGetConvergenceStatus_HIOPSPARSEGPU(OPFLOW opflow, + PetscBool *status) { + OPFLOWSolver_HIOPSPARSEGPU hiop = (OPFLOWSolver_HIOPSPARSEGPU)opflow->solver; + + PetscFunctionBegin; + if (hiop->status < 3) + *status = PETSC_TRUE; /* See hiopInterface.hpp. The first three denote + convergence */ + else + *status = PETSC_FALSE; + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverGetObjective_HIOPSPARSEGPU(OPFLOW opflow, + PetscReal *obj) { + PetscFunctionBegin; + *obj = opflow->obj; + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverGetSolution_HIOPSPARSEGPU(OPFLOW opflow, Vec *X) { + PetscFunctionBegin; + *X = opflow->X; + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverGetConstraints_HIOPSPARSEGPU(OPFLOW opflow, Vec *G) { + PetscFunctionBegin; + *G = opflow->G; + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverGetConstraintMultipliers_HIOPSPARSEGPU(OPFLOW opflow, + Vec *Lambda) { + PetscFunctionBegin; + *Lambda = opflow->Lambda; + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverDestroy_HIOPSPARSEGPU(OPFLOW opflow) { + PetscErrorCode ierr; + OPFLOWSolver_HIOPSPARSEGPU hiop = (OPFLOWSolver_HIOPSPARSEGPU)opflow->solver; + + PetscFunctionBegin; + + delete hiop->solver; + delete hiop->sp; + delete hiop->nlp; + + ierr = PetscFree(hiop); + CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +PetscErrorCode OPFLOWSolverCreate_HIOPSPARSEGPU(OPFLOW opflow) { + PetscErrorCode ierr; + OPFLOWSolver_HIOPSPARSEGPU hiop; + + PetscFunctionBegin; + ierr = PetscCalloc1(1, &hiop); + CHKERRQ(ierr); + + opflow->solver = hiop; + + opflow->solverops.setup = OPFLOWSolverSetUp_HIOPSPARSEGPU; + opflow->solverops.solve = OPFLOWSolverSolve_HIOPSPARSEGPU; + opflow->solverops.destroy = OPFLOWSolverDestroy_HIOPSPARSEGPU; + opflow->solverops.getobjective = OPFLOWSolverGetObjective_HIOPSPARSEGPU; + opflow->solverops.getconvergencestatus = + OPFLOWSolverGetConvergenceStatus_HIOPSPARSEGPU; + opflow->solverops.getsolution = OPFLOWSolverGetSolution_HIOPSPARSEGPU; + opflow->solverops.getconstraints = OPFLOWSolverGetConstraints_HIOPSPARSEGPU; + opflow->solverops.getconstraintmultipliers = + OPFLOWSolverGetConstraintMultipliers_HIOPSPARSEGPU; + + PetscFunctionReturn(0); +} + +#endif +#endif diff --git a/src/opflow/solver/hiop/opflow_hiopsparsegpu.hpp b/src/opflow/solver/hiop/opflow_hiopsparsegpu.hpp new file mode 100644 index 00000000..de81eb4c --- /dev/null +++ b/src/opflow/solver/hiop/opflow_hiopsparsegpu.hpp @@ -0,0 +1,120 @@ +#include + +#if defined(EXAGO_ENABLE_HIOP) +#if defined(EXAGO_ENABLE_HIOP_SPARSE) + +#ifndef OPFLOWHIOPSPARSEGPU_HPP +#define OPFLOWHIOPSPARSEGPU_HPP +#include +#include +#include +#include +#include + +#if defined(EXAGO_ENABLE_IPOPT) +#include // ipopt adapter from hiop +#include +#endif + +#include +#include +#include +#include + +#ifdef HIOP_USE_MAGMA +#include "magma_v2.h" +#endif + +class OPFLOWHIOPSPARSEGPUInterface : public hiop::hiopInterfaceSparse { +public: + OPFLOWHIOPSPARSEGPUInterface(OPFLOW); + + ~OPFLOWHIOPSPARSEGPUInterface() {} + + bool get_prob_sizes(hiop::size_type &n, hiop::size_type &m); + + bool get_vars_info(const hiop::size_type &n, double *xlow, double *xupp, + NonlinearityType *type); + + bool get_cons_info(const hiop::size_type &m, double *clow, double *cupp, + NonlinearityType *type); + + bool get_sparse_blocks_info(hiop::size_type &nx, + hiop::size_type &nnz_sparse_Jaceq, + hiop::size_type &nnz_sparse_Jacineq, + hiop::size_type &nnz_sparse_Hess_Lagr); + + bool eval_f(const hiop::size_type &n, const double *x, bool new_x, + double &obj_value); + + bool eval_cons(const hiop::size_type &n, const hiop::size_type &m, + const hiop::size_type &num_cons, + const hiop::size_type *idx_cons, const double *x, bool new_x, + double *cons); + + bool eval_cons(const hiop::size_type &n, const hiop::size_type &m, + const double *x, bool new_x, double *cons); + + bool eval_grad_f(const hiop::size_type &n, const double *x, bool new_x, + double *gradf); + + bool eval_Jac_cons(const hiop::size_type &n, const hiop::size_type &m, + const hiop::size_type &num_cons, + const hiop::index_type *idx_cons, const double *x, + bool new_x, const hiop::size_type &nnzJacS, + hiop::index_type *iJacS, hiop::index_type *jJacS, + double *MJacS); + + bool eval_Jac_cons(const hiop::size_type &n, const hiop::size_type &m, + const double *x, bool new_x, + const hiop::size_type &nnzJacS, hiop::index_type *iJacS, + hiop::index_type *jJacS, double *MJacS); + + bool get_starting_point(const hiop::size_type &n, double *x0); + bool eval_Hess_Lagr(const hiop::size_type &n, const hiop::size_type &m, + const double *x, bool new_x, const double &obj_factor, + const double *lambda, bool new_lambda, + const hiop::size_type &nnzHSS, hiop::index_type *iHSS, + hiop::index_type *jHSS, double *MHSS); + + void solution_callback(hiop::hiopSolveStatus status, int n, const double *x, + const double *z_L, const double *z_U, int m, + const double *g, const double *lambda, + double obj_value); + + bool iterate_callback(int iter, double obj_value, double logbar_obj_value, + int n, const double *x, const double *z_L, + const double *z_U, int m_ineq, const double *s, int m, + const double *g, const double *lambda, double inf_pr, + double inf_du, double onenorm_pr_, double mu, + double alpha_du, double alpha_pr, int ls_trials); + + bool get_MPI_comm(MPI_Comm &comm_out) { + comm_out = MPI_COMM_SELF; + return true; + } + +private: + OPFLOW opflow; +}; + +typedef struct _p_OPFLOWSolver_HIOPSPARSEGPU *OPFLOWSolver_HIOPSPARSEGPU; + +struct _p_OPFLOWSolver_HIOPSPARSEGPU { + + OPFLOWHIOPSPARSEGPUInterface *nlp; + hiop::hiopSolveStatus status; + hiop::hiopNlpSparse *sp; + hiop::hiopAlgFilterIPMNewton *solver; + +#if defined(EXAGO_ENABLE_IPOPT) + // Ipopt Adapter structs + SmartPtr ipoptTNLP; + SmartPtr ipoptApp; + PetscBool ipopt_debug; +#endif +}; + +#endif +#endif +#endif diff --git a/tests/functionality/opflow/CMakeLists.txt b/tests/functionality/opflow/CMakeLists.txt index 10ea9d79..aefe3dac 100644 --- a/tests/functionality/opflow/CMakeLists.txt +++ b/tests/functionality/opflow/CMakeLists.txt @@ -64,6 +64,26 @@ if(EXAGO_INSTALL_TESTS) ) endif() + if(EXAGO_ENABLE_GPU) + exago_add_test( + NAME + FUNCTIONALITY_TEST_OPFLOW_RAJAHIOP_SPARSE_GPU_TOML_TESTSUITE + COMMAND + ${RUNCMD} + $ + ${CMAKE_CURRENT_SOURCE_DIR}/hiop_pbpolrajahiopsparse_gpu.toml + DEPENDS + HIOP + RAJA + ) + if(EXAGO_ENABLE_HIOP AND EXAGO_ENABLE_RAJA) + set_tests_properties( + FUNCTIONALITY_TEST_OPFLOW_RAJAHIOP_SPARSE_GPU_TOML_TESTSUITE + PROPERTIES LABELS "incline-skip;deception-skip;newell-skip" + ) + endif() + endif(EXAGO_ENABLE_GPU) + # Note: All cartesian and/or current balance models are way behind. These # models should not be used and their tests should be disabled. if(EXAGO_FEATURE_CARTESIAN_ACTIVE) diff --git a/tests/functionality/opflow/hiop_pbpolrajahiopsparse_gpu.toml b/tests/functionality/opflow/hiop_pbpolrajahiopsparse_gpu.toml new file mode 100644 index 00000000..c509b0c5 --- /dev/null +++ b/tests/functionality/opflow/hiop_pbpolrajahiopsparse_gpu.toml @@ -0,0 +1,47 @@ +testsuite_name = "HiOp OPFLOW Functionality" + +[presets] +solver = 'HIOPSPARSEGPU' +model = 'PBPOLRAJAHIOPSPARSE' +hiop_compute_mode = 'GPU' +hiop_mem_space = 'DEVICE' +tolerance = 1e-6 +warning_tolerance = 0.01 +gen_bus_voltage_type = 'VARIABLE_WITHIN_BOUNDS' +has_gen_set_point = false +is_loadloss_active = false +is_powerimb_active = false +use_agc = false +load_loss_penalty = 1000.0 +power_imbalance_penalty = 1000.0 +initialization_type = 'MIDPOINT' +iter_range = 0 + +[[testcase]] +network = 'datafiles/case9/case9mod.m' +description = 'datafiles/case9/case9mod.m base case' +is_loadloss_active = false +is_powerimb_active = false +num_iters = 15 #usually 10 or 11 +iter_range = 1 +obj_value = 4.1444507e+03 + +[[testcase]] +network = 'datafiles/case9/case9mod.m' +description = "datafiles/case9/case9mod.m initialization DCOPF" +is_loadloss_active = false +is_powerimb_active = false +initialization_type = "DCOPF" +num_iters = 12 +iter_range = 1 +obj_value = 4.1444507e+03 + +[[testcase]] +network = 'datafiles/case_ACTIVSg200.m' +description = "datafiles/case_ACTIVSg200.m initialization DCOPF" +is_loadloss_active = false +is_powerimb_active = false +initialization_type = "DCOPF" +num_iters = 38 #usually 16 or 31 +iter_range = 15 +obj_value = 2.7552967e+04