Skip to content

Commit d92d4ff

Browse files
lochelarun3688
lochel
and
arun3688
authored
Use fmi2_import_get_directional_derivative to construct Jacobian matrix (#1122)
Co-authored-by: arun3688 <rain100falls@gmail.com>
1 parent 495b3a7 commit d92d4ff

File tree

15 files changed

+211
-125
lines changed

15 files changed

+211
-125
lines changed

‎src/OMSimulatorLib/AlgLoop.cpp‎

Lines changed: 119 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -50,9 +50,10 @@ inline bool checkFlag(int flag, std::string functionName)
5050
{
5151
if (flag < 0)
5252
{
53-
logError("SUNDIALS_ERROR: " + functionName + "() failed with flag = " + std::to_string(flag));
53+
logError("SUNDIALS_ERROR: " + functionName + " failed with flag = " + std::to_string(flag));
5454
return false;
5555
}
56+
logDebug("SUNDIALS_INFO: " + functionName + " failed with flag = " + std::to_string(flag));
5657
return true;
5758
}
5859

@@ -63,25 +64,25 @@ inline bool checkFlag(int flag, std::string functionName)
6364
* @param module Name of the module reporting the error.
6465
* @param function Name of the function in which the error occurred.
6566
* @param msg Error Message.
66-
* @param userData Pointer to user data. Unused.
67+
* @param user_data Pointer to user data. Unused.
6768
*/
6869
void oms::KinsolSolver::sundialsErrorHandlerFunction(int errorCode, const char *module,
6970
const char *function, char *msg,
70-
void *userData)
71+
void *user_data)
7172
{
7273
KINSOL_USER_DATA* kinsolUserData;
7374
std::string systNum = "unknown";
7475
std::string mod = module;
7576
std::string func = function;
7677

77-
if (userData != NULL) {
78-
kinsolUserData = (KINSOL_USER_DATA *)userData;
78+
if (user_data != NULL)
79+
{
80+
kinsolUserData = (KINSOL_USER_DATA *)user_data;
7981
systNum = std::to_string(kinsolUserData->algLoopNumber);
8082
}
8183

8284
logError("SUNDIALS_ERROR: [system] " + systNum + " [module] " + mod + " | [function] " + func
83-
+" | [error_code] " + std::to_string(errorCode));
84-
logError(msg);
85+
+ " | [error_code] " + std::to_string(errorCode) + "\n" + std::string(msg));
8586
}
8687

8788
/**
@@ -92,39 +93,106 @@ void oms::KinsolSolver::sundialsErrorHandlerFunction(int errorCode, const char *
9293
* @param module Name of the module reporting the information.
9394
* @param function Name of the function reporting the information.
9495
* @param msg Message.
95-
* @param userData Pointer to user data. Unused.
96+
* @param user_data Pointer to user data. Unused.
9697
*/
9798
void oms::KinsolSolver::sundialsInfoHandlerFunction(const char *module, const char *function,
98-
char *msg, void *userData)
99+
char *msg, void *user_data)
99100
{
100101
KINSOL_USER_DATA* kinsolUserData;
101102
std::string systNum = "unknown";
102103
std::string mod = module;
103104
std::string func = function;
104105

105-
if (userData != NULL) {
106-
kinsolUserData = (KINSOL_USER_DATA *)userData;
106+
if (user_data != NULL) {
107+
kinsolUserData = (KINSOL_USER_DATA *)user_data;
107108
systNum = std::to_string(kinsolUserData->algLoopNumber);
108109
}
109110

110-
logDebug("SUNDIALS_INFO: [system] " + systNum + " [module] " + mod + " | [function] " + func);
111-
logDebug(msg);
111+
logDebug("SUNDIALS_INFO: [system] " + systNum + " [module] " + mod + " | [function] " + func + "\n" + std::string(msg));
112+
}
113+
114+
/**
115+
* @brief Jacobian function for KINSOL
116+
*
117+
* @param u is the current (unscaled) iterate
118+
* @param fu is the current value of the vector F(u)
119+
* @param J is the output approximate Jacobian matrix, J = ∂F/∂u,
120+
* of type SUNMatrix
121+
* @param data is a pointer to user data, the same as the user data
122+
* parameter passed to KINSetUserData
123+
* @param tmp1 pointer to memory allocated for variables of type
124+
* N_Vector which can be used by the KINJacFn function as
125+
* temporary storage or work space
126+
* @param tmp2 pointer to memory allocated for variables of type
127+
* N_Vector which can be used by the KINJacFn function as
128+
* temporary storage or work space
129+
* @return int A function of type KINLsJacFn should return 0 if
130+
* successful or a non-zero value otherwise
131+
*/
132+
int oms::KinsolSolver::nlsKinsolJac(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
133+
{
134+
KINSOL_USER_DATA* kinsoluserData = (KINSOL_USER_DATA*) user_data;
135+
System* syst = kinsoluserData->syst;
136+
AlgLoop* algLoop = syst->getAlgLoop(kinsoluserData->algLoopNumber);
137+
DirectedGraph* graph = kinsoluserData->graph;
138+
const oms_ssc_t SCC = algLoop->getSCC();
139+
const int size = SCC.size();
140+
141+
oms_status_enu_t status;
142+
double *u_data = NV_DATA_S(u);
143+
144+
for(int j=0; j<size; j++)
145+
{
146+
const int output = SCC[j].first;
147+
const oms::ComRef& crefOutput = graph->getNodes()[output].getName();
148+
const int input = SCC[j].second;
149+
const oms::ComRef& crefInput = graph->getNodes()[input].getName();
150+
// res[j] = output - input
151+
// std::cout << "res[" << j << "] = " << std::string(crefOutput) << " - " << std::string(crefInput) << std::endl;
152+
153+
for(int i=0; i<size; i++)
154+
{
155+
double der = 0.0;
156+
const int input_col = SCC[i].second;
157+
oms::ComRef crefInput_col = graph->getNodes()[input_col].getName();
158+
oms::ComRef front = crefInput_col.pop_front();
159+
160+
if (front == crefOutput.front())
161+
if (oms_status_ok != syst->getDirectionalDerivative(crefOutput, crefInput_col, der))
162+
return logError("not recoverable error");
163+
if (input_col == input)
164+
der -= 1.0;
165+
166+
SM_ELEMENT_D(J, j, i) = der;
167+
}
168+
}
169+
170+
return 0;
112171
}
113172

114173
/**
115174
* @brief Residual function for KINSOL
116175
*
117-
* @param uu Input value
118-
* @param fval Contains residual at output
119-
* @param userData Pointer to user data to access System, AlgLoop, Directed Graph and CSS as well as getReal and setReal
120-
* @return int Return 0 on success, 1 if an recoverable error occured, -1 if a fatal error occured.
176+
* This function computes F(u) (or G(u) for fixed-point iteration and
177+
* Anderson acceleration) for a given value of the vector u.
178+
*
179+
* @param u is the current value of the variable vector, u
180+
* @param fval is the output vector F(u)
181+
* @param user_data a is a pointer to user data, the pointer user
182+
* data passed to KINSetUserData
183+
* @return int A KINSysFn function should return 0 if
184+
* successful, a positive value if a recoverable
185+
* error occurred (in which case kinsol will attempt
186+
* to correct), or a negative value if it failed
187+
* unrecoverably (in which case the solution process
188+
* is halted and KIN SYSFUNC FAIL is returned).
121189
*/
122-
int oms::KinsolSolver::nlsKinsolResiduals(N_Vector uu, N_Vector fval, void *userData)
190+
int oms::KinsolSolver::nlsKinsolResiduals(N_Vector u, N_Vector fval, void *user_data)
123191
{
124-
double *uu_data = NV_DATA_S(uu);
192+
double *u_data = NV_DATA_S(u);
125193
double *fval_data = NV_DATA_S(fval);
126194

127-
KINSOL_USER_DATA* kinsoluserData = (KINSOL_USER_DATA*) userData;
195+
KINSOL_USER_DATA* kinsoluserData = (KINSOL_USER_DATA*) user_data;
128196
System* syst = kinsoluserData->syst;
129197
AlgLoop* algLoop = syst->getAlgLoop(kinsoluserData->algLoopNumber);
130198
DirectedGraph* graph = kinsoluserData->graph;
@@ -142,13 +210,13 @@ int oms::KinsolSolver::nlsKinsolResiduals(N_Vector uu, N_Vector fval, void *user
142210
ss << "inputs:" << std::endl;
143211
}
144212

145-
// Set values from uu
213+
// Set values from u
146214
for (int i=0; i<size; ++i)
147215
{
148-
int output = SCC[i].second;
149-
status = syst->setReal(graph->getNodes()[output].getName(), uu_data[i]);
216+
int input = SCC[i].second;
217+
status = syst->setReal(graph->getNodes()[input].getName(), u_data[i]);
150218
if (Flags::DumpAlgLoops())
151-
ss << " " << graph->getNodes()[output].getName().c_str() << ": " << uu_data[i] << std::endl;
219+
ss << " " << graph->getNodes()[input].getName().c_str() << ": " << u_data[i] << std::endl;
152220
if (status == oms_status_discard || status == oms_status_error || status == oms_status_warning)
153221
{
154222
logInfo("iteration " + std::to_string(kinsoluserData->iteration) + ": recoverable error (1)");
@@ -181,7 +249,7 @@ int oms::KinsolSolver::nlsKinsolResiduals(N_Vector uu, N_Vector fval, void *user
181249
logInfo("iteration " + std::to_string(kinsoluserData->iteration) + ": not recoverable error (2)");
182250
return -1 /* not recoverable error */;
183251
}
184-
fval_data[i] = fval_data[i] - uu_data[i];
252+
fval_data[i] = fval_data[i] - u_data[i];
185253
}
186254

187255
if (Flags::DumpAlgLoops())
@@ -212,7 +280,7 @@ oms::KinsolSolver::~KinsolSolver()
212280
SUNMatDestroy(this->J);
213281
N_VDestroy_Serial(this->y);
214282

215-
delete((KINSOL_USER_DATA*)(this->userData));
283+
delete((KINSOL_USER_DATA*)(this->user_data));
216284
}
217285

218286
/**
@@ -223,7 +291,7 @@ oms::KinsolSolver::~KinsolSolver()
223291
* @param absoluteTolerance Tolerance used for solving the loop
224292
* @return oms::KinsolSolver* Retruns pointer to KinsolSolver object
225293
*/
226-
oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance)
294+
oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, const bool useDirectionalDerivative)
227295
{
228296
int flag;
229297
int printLevel;
@@ -250,24 +318,27 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
250318
}
251319

252320
/* Set user data given to KINSOL */
253-
kinsolSolver->userData = new KINSOL_USER_DATA{/*.syst=*/NULL, /*.graph=*/NULL, /*.algLoopNumber=*/algLoopNum, /*.iteration=*/0};
254-
flag = KINSetUserData(kinsolSolver->kinsolMemory, kinsolSolver->userData);
321+
kinsolSolver->user_data = new KINSOL_USER_DATA{/*.syst=*/NULL, /*.graph=*/NULL, /*.algLoopNumber=*/algLoopNum, /*.iteration=*/0};
322+
flag = KINSetUserData(kinsolSolver->kinsolMemory, kinsolSolver->user_data);
255323
if (!checkFlag(flag, "KINSetUserData")) return NULL;
256324

257325
/* Set error handler and print level */
258-
if (Log::DebugEnabled()) {
326+
if (Log::DebugEnabled())
327+
{
259328
logDebug("SUNDIALS KINSOL: Set print level to maximum.");
260329
printLevel = 3;
261-
} else {
330+
}
331+
else
332+
{
262333
printLevel = 0;
263334
}
264335
flag = KINSetPrintLevel(kinsolSolver->kinsolMemory, printLevel);
265336
if (!checkFlag(flag, "KINSetPrintLevel")) return NULL;
266337

267-
flag = KINSetErrHandlerFn(kinsolSolver->kinsolMemory, sundialsErrorHandlerFunction, kinsolSolver->userData);
338+
flag = KINSetErrHandlerFn(kinsolSolver->kinsolMemory, sundialsErrorHandlerFunction, kinsolSolver->user_data);
268339
if (!checkFlag(flag, "KINSetErrHandlerFn")) return NULL;
269340

270-
flag = KINSetInfoHandlerFn(kinsolSolver->kinsolMemory, sundialsInfoHandlerFunction, kinsolSolver->userData);
341+
flag = KINSetInfoHandlerFn(kinsolSolver->kinsolMemory, sundialsInfoHandlerFunction, kinsolSolver->user_data);
271342
if (!checkFlag(flag, "KINSetInfoHandlerFn")) return NULL;
272343

273344
/* Initialize KINSOL object */
@@ -285,7 +356,10 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
285356
if (!checkFlag(flag, "KINSetLinearSolver")) return NULL;
286357

287358
/* Set Jacobian for linear solver */
288-
flag = KINSetJacFn(kinsolSolver->kinsolMemory, NULL); /* Use KINSOLs internal difference quotient function for Jacobian approximation */
359+
if (useDirectionalDerivative && Flags::DirectionalDerivatives())
360+
flag = KINSetJacFn(kinsolSolver->kinsolMemory, nlsKinsolJac); /* Use symbolic Jacobian */
361+
else
362+
flag = KINSetJacFn(kinsolSolver->kinsolMemory, NULL); /* Use KINSOLs internal difference quotient function for Jacobian approximation */
289363
if (!checkFlag(flag, "KINSetJacFn")) return NULL;
290364

291365
/* Set function-norm stopping tolerance */
@@ -328,7 +402,7 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
328402
oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& graph)
329403
{
330404
/* Update user data */
331-
KINSOL_USER_DATA* kinsolUserData = (KINSOL_USER_DATA*) userData;
405+
KINSOL_USER_DATA* kinsolUserData = (KINSOL_USER_DATA*) user_data;
332406
kinsolUserData->syst = &syst;
333407
kinsolUserData->graph = &graph;
334408
kinsolUserData->iteration = 0;
@@ -364,24 +438,16 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
364438
flag = KINSol(kinsolMemory, /* KINSol memory block */
365439
initialGuess, /* initial guess on input; solution vector */
366440
KIN_NONE, /* global strategy choice: Basic newton iteration */
367-
uScale, /* scaling vector, for the variable uu */
441+
uScale, /* scaling vector, for the variable u */
368442
fScale); /* scaling vector for function values fval */
369-
if (flag < 0)
370-
{
371-
logError("SUNDIALS_ERROR: KINSol() failed with flag = " + std::to_string(flag));
372-
return oms_status_error;
373-
}
374-
else
375-
{
376-
logDebug("SUNDIALS_INFO: KINSol() succeded with flag = " + std::to_string(flag));
377-
}
443+
if (!checkFlag(flag, "KINSol")) return oms_status_error;
378444

379445
/* Check solution */
380-
flag = nlsKinsolResiduals(initialGuess, fTmp, userData);
446+
flag = nlsKinsolResiduals(initialGuess, fTmp, user_data);
381447
fNormValue = N_VWL2Norm(fTmp, fTmp);
382448
if ( fNormValue > fnormtol )
383449
{
384-
logWarning("Solution of algebraic loop " + std::to_string(((KINSOL_USER_DATA *)userData)->algLoopNumber) + "not within precission given by fnormtol: " + std::to_string(fnormtol));
450+
logWarning("Solution of algebraic loop " + std::to_string(((KINSOL_USER_DATA *)user_data)->algLoopNumber) + "not within precission given by fnormtol: " + std::to_string(fnormtol));
385451
logDebug("2-norm of residual of solution: " + std::to_string(fNormValue));
386452
return oms_status_warning;
387453
}
@@ -400,7 +466,7 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
400466
* @param scc Strong Connected Componten, a vector of connected
401467
* @param number
402468
*/
403-
oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t scc, const int number): absoluteTolerance(absTol), SCC(scc), systNumber(number)
469+
oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t scc, const int number, const bool useDirectionalDerivative): absoluteTolerance(absTol), SCC(scc), systNumber(number)
404470
{
405471
switch (method)
406472
{
@@ -415,7 +481,7 @@ oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t scc,
415481

416482
if (method == oms_alg_solver_kinsol)
417483
{
418-
kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.size(), absoluteTolerance);
484+
kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.size(), absoluteTolerance, useDirectionalDerivative);
419485
if (kinsolData==NULL)
420486
{
421487
logError("NewKinsolSolver() failed. Aborting!");
@@ -574,18 +640,20 @@ std::string oms::AlgLoop::getAlgSolverName(void)
574640
std::string oms::AlgLoop::dumpLoopVars(DirectedGraph& graph)
575641
{
576642
const int size = SCC.size();
577-
std::string varNames = "\t";
643+
std::string varNames = "";
578644
int output;
579645

580646
for (int i=0; i<size-1; ++i)
581647
{
648+
varNames.append(" ");
582649
output = SCC[i].first;
583650
varNames.append(graph.getNodes()[output].getName().c_str());
584651
varNames.append(" -> ");
585652
output = SCC[i].second;
586653
varNames.append(graph.getNodes()[output].getName().c_str());
587-
varNames.append("\n\t");
654+
varNames.append("\n");
588655
}
656+
varNames.append(" ");
589657
output = SCC[size-1].first;
590658
varNames.append(graph.getNodes()[output].getName().c_str());
591659
varNames.append(" -> ");

‎src/OMSimulatorLib/AlgLoop.h‎

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -55,12 +55,12 @@ namespace oms
5555

5656
class KinsolSolver
5757
{
58-
public:
58+
public:
5959
~KinsolSolver();
60-
static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance);
60+
static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, const bool useDirectionalDerivative);
6161
oms_status_enu_t kinsolSolve(System& syst, DirectedGraph& graph);
6262

63-
private:
63+
private:
6464
/* tolerances */
6565
double fnormtol; /* function tolerance */
6666

@@ -72,7 +72,7 @@ namespace oms
7272

7373
/* kinsol internal data */
7474
void* kinsolMemory;
75-
void* userData;
75+
void* user_data;
7676
int size;
7777

7878
/* linear solver data */
@@ -81,22 +81,23 @@ namespace oms
8181
SUNMatrix J; /* (Non-)Sparse matrix template for cloning matrices needed within linear solver */
8282

8383
/* member function */
84-
static int nlsKinsolResiduals(N_Vector uu, N_Vector fval, void *userData);
85-
static void sundialsErrorHandlerFunction(int error_code, const char *module, const char *function, char *msg, void *userData);
86-
static void sundialsInfoHandlerFunction(const char *module, const char *function, char *msg, void *userData);
84+
static int nlsKinsolJac(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2);
85+
static int nlsKinsolResiduals(N_Vector u, N_Vector fval, void *user_data);
86+
static void sundialsErrorHandlerFunction(int error_code, const char *module, const char *function, char *msg, void *user_data);
87+
static void sundialsInfoHandlerFunction(const char *module, const char *function, char *msg, void *user_data);
8788
};
8889

8990
class AlgLoop
9091
{
91-
public:
92-
AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t SCC, const int systNumber);
92+
public:
93+
AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t SCC, const int systNumber, const bool useDirectionalDerivative);
9394

9495
oms_ssc_t getSCC() {return SCC;}
9596
oms_status_enu_t solveAlgLoop(System& syst, DirectedGraph& graph);
9697
std::string getAlgSolverName();
9798
std::string dumpLoopVars(DirectedGraph& graph);
9899

99-
private:
100+
private:
100101
oms_alg_solver_enu_t algSolverMethod;
101102
oms_status_enu_t fixPointIteration(System& syst, DirectedGraph& graph);
102103

@@ -109,4 +110,4 @@ namespace oms
109110
};
110111
}
111112

112-
#endif // _OMS_ALGLOOP_H_
113+
#endif

0 commit comments

Comments
 (0)