OMCompiler/SimulationRuntime/c/simulation/eval_dep.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * This file belongs to the OpenModelica Run-Time System | ||
| 3 | * | ||
| 4 | * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings | ||
| 5 | * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights | ||
| 6 | * reserved. | ||
| 7 | * | ||
| 8 | * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE | ||
| 9 | * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY | ||
| 10 | * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S | ||
| 11 | * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL | ||
| 12 | * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. | ||
| 13 | * | ||
| 14 | * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License | ||
| 15 | * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs: | ||
| 16 | * http://www.openmodelica.org or https://github.com/OpenModelica/ or | ||
| 17 | * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU | ||
| 18 | * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW | ||
| 19 | * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause. | ||
| 20 | * | ||
| 21 | * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 22 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY | ||
| 23 | * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF | ||
| 24 | * OSMC-PL. | ||
| 25 | * | ||
| 26 | */ | ||
| 27 | |||
| 28 | #include "eval_dep.h" | ||
| 29 | #include "../util/omc_error.h" | ||
| 30 | #include "../util/uthash.h" | ||
| 31 | #include "../simulation_data.h" | ||
| 32 | #include "simulation_info_json.h" | ||
| 33 | |||
| 34 | |||
| 35 | /** | ||
| 36 | * @brief allocate new empty DAG with no edges | ||
| 37 | * | ||
| 38 | * @param nVars number of variables in the system | ||
| 39 | * @param nEqns number of equations in the system | ||
| 40 | * @return EVAL_DAG* | ||
| 41 | */ | ||
| 42 | ✗ | EVAL_DAG* allocEvalDAG(size_t nVars, size_t nEqns) | |
| 43 | { | ||
| 44 | ✗ | EVAL_DAG* dag = (EVAL_DAG*) malloc(sizeof(EVAL_DAG)); | |
| 45 | ✗ | dag->nVars = nVars; | |
| 46 | ✗ | dag->nEqns = nEqns; | |
| 47 | ✗ | if (nVars != 0) { | |
| 48 | ✗ | dag->mapVarToEqNode = (size_t*) malloc(nVars*sizeof(size_t)); | |
| 49 | } else { | ||
| 50 | ✗ | dag->mapVarToEqNode = NULL; | |
| 51 | } | ||
| 52 | ✗ | if (nEqns != 0 ) { | |
| 53 | ✗ | dag->nEqDep = (size_t*) calloc(nEqns, sizeof(size_t)); | |
| 54 | ✗ | dag->eqDep = (size_t**) calloc(nEqns, sizeof(size_t*)); | |
| 55 | ✗ | dag->select = (int*) calloc(nEqns, sizeof(int)); | |
| 56 | } else { | ||
| 57 | ✗ | dag->nEqDep = NULL; | |
| 58 | ✗ | dag->eqDep = NULL; | |
| 59 | ✗ | dag->select = NULL; | |
| 60 | } | ||
| 61 | |||
| 62 | /* | ||
| 63 | workaround, some JACOBIAN_TMP_VAR variables seem to miss an equation in which | ||
| 64 | they are solved. This is probably because the derivative is zero and the | ||
| 65 | corresponding equation was removed from the system. | ||
| 66 | |||
| 67 | So we use the default -1 to mean there is no equation to evaluate. | ||
| 68 | |||
| 69 | A better solution would be to remove the variable as well and propagate the | ||
| 70 | zero symbolically. | ||
| 71 | */ | ||
| 72 | ✗ | for (size_t i = 0; i < nVars; ++i) { | |
| 73 | ✗ | dag->mapVarToEqNode[i] = (size_t)(-1); | |
| 74 | } | ||
| 75 | |||
| 76 | ✗ | return dag; | |
| 77 | } | ||
| 78 | |||
| 79 | /** | ||
| 80 | * @brief Free DAG | ||
| 81 | * | ||
| 82 | * @param dag | ||
| 83 | */ | ||
| 84 | 7 | void freeEvalDAG(EVAL_DAG* dag) | |
| 85 | { | ||
| 86 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
|
7 | if (dag) { |
| 87 | ✗ | free(dag->select); | |
| 88 | /* only free once since eqDep was allocated in one chunk */ | ||
| 89 | ✗ | if (dag->eqDep) free(dag->eqDep[0]); | |
| 90 | ✗ | free(dag->eqDep); | |
| 91 | ✗ | free(dag->nEqDep); | |
| 92 | ✗ | free(dag->mapVarToEqNode); | |
| 93 | ✗ | free(dag); | |
| 94 | } | ||
| 95 | 7 | } | |
| 96 | |||
| 97 | /** | ||
| 98 | * @brief Create empty selection | ||
| 99 | * | ||
| 100 | * @param dag | ||
| 101 | * @return EVAL_SELECTION* | ||
| 102 | */ | ||
| 103 | ✗ | EVAL_SELECTION* allocEvalSelection(EVAL_DAG* dag) | |
| 104 | { | ||
| 105 | ✗ | assertStreamPrint(NULL, dag, "No DAG was given."); | |
| 106 | |||
| 107 | ✗ | EVAL_SELECTION* selection = (EVAL_SELECTION*) malloc(sizeof(EVAL_SELECTION)); | |
| 108 | ✗ | selection->n = 0; | |
| 109 | ✗ | selection->idx = (size_t*) malloc(dag->nEqns * sizeof(size_t)); | |
| 110 | ✗ | selection->dag = dag; | |
| 111 | |||
| 112 | ✗ | return selection; | |
| 113 | } | ||
| 114 | |||
| 115 | /** | ||
| 116 | * @brief Free selection | ||
| 117 | * | ||
| 118 | * except the DAG, it may be used by other selections. | ||
| 119 | * | ||
| 120 | * @param selection | ||
| 121 | */ | ||
| 122 | 6 | void freeEvalSelection(EVAL_SELECTION* selection) | |
| 123 | { | ||
| 124 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (selection) { |
| 125 | ✗ | free(selection->idx); | |
| 126 | ✗ | free(selection); | |
| 127 | } | ||
| 128 | 6 | } | |
| 129 | |||
| 130 | /** | ||
| 131 | * @brief clear selection | ||
| 132 | * | ||
| 133 | * call this first, then select manually, then activateEvalDependencies | ||
| 134 | * | ||
| 135 | * @param selection | ||
| 136 | */ | ||
| 137 | ✗ | void clearEvalSelection(EVAL_SELECTION* selection) | |
| 138 | { | ||
| 139 | ✗ | assertStreamPrint(NULL, selection, "selection is NULL."); | |
| 140 | |||
| 141 | /* clear work array */ | ||
| 142 | ✗ | for (size_t i = 0; i < selection->dag->nEqns; ++i) { | |
| 143 | ✗ | selection->dag->select[i] = 0 /* FALSE */; | |
| 144 | } | ||
| 145 | |||
| 146 | /* set selected equations to zero just to be safe */ | ||
| 147 | ✗ | selection->n = 0; | |
| 148 | // don't clear idx as that would be O(n) work, | ||
| 149 | // it will be overwritten by activateEvalDependencies | ||
| 150 | ✗ | } | |
| 151 | |||
| 152 | /** | ||
| 153 | * @brief Set dependencies based on already selected subset. | ||
| 154 | * | ||
| 155 | * Since we have a DAG we can just go backwards and look at direct dependency | ||
| 156 | * only, because indirect dependencies will be handled when we get to the direct | ||
| 157 | * dependency node which sets its direct dependency and so on... | ||
| 158 | */ | ||
| 159 | ✗ | void activateEvalDependencies(EVAL_SELECTION* selection) | |
| 160 | { | ||
| 161 | ✗ | assertStreamPrint(NULL, selection, "selection is NULL."); | |
| 162 | |||
| 163 | ✗ | EVAL_DAG* dag = selection->dag; | |
| 164 | /* select dependencies backwards */ | ||
| 165 | ✗ | for (size_t i = dag->nEqns-1; i+1 > 0 /* careful: count down with unsigned */; --i) { | |
| 166 | ✗ | if (dag->select[i]) { | |
| 167 | ✗ | for (size_t j = 0; j < dag->nEqDep[i]; ++j) { | |
| 168 | ✗ | size_t dep = dag->eqDep[i][j]; | |
| 169 | |||
| 170 | /* | ||
| 171 | workaround, some JACOBIAN_TMP_VAR variables seem to miss an equation in | ||
| 172 | which they are solved. This is probably because the derivative is zero | ||
| 173 | and the corresponding equation was removed from the system. | ||
| 174 | |||
| 175 | So we use the default -1 to mean there is no equation to evaluate. | ||
| 176 | |||
| 177 | A better solution would be to remove the variable as well and propagate | ||
| 178 | the zero symbolically. | ||
| 179 | */ | ||
| 180 | ✗ | if (dep == (size_t)(-1)) continue; | |
| 181 | |||
| 182 | ✗ | dag->select[dep] = !0 /* TRUE */; | |
| 183 | } | ||
| 184 | } | ||
| 185 | } | ||
| 186 | |||
| 187 | /* get the indices in correct order */ | ||
| 188 | ✗ | selection->n = 0; | |
| 189 | ✗ | for (size_t i = 0; i < dag->nEqns; ++i) { | |
| 190 | ✗ | if (dag->select[i]) selection->idx[selection->n++] = i; | |
| 191 | } | ||
| 192 | ✗ | } | |
| 193 | |||
| 194 | /* * * * * * * * * * * * * * * | ||
| 195 | * OpenModelica specific stuff | ||
| 196 | * * * * * * * * * * * * * * */ | ||
| 197 | |||
| 198 | typedef struct hash_varName_index | ||
| 199 | { | ||
| 200 | const char *id; // variable name | ||
| 201 | size_t varIndex; // variable index | ||
| 202 | size_t eqIndex; // equation index | ||
| 203 | UT_hash_handle hh; | ||
| 204 | } hash_varName_index; | ||
| 205 | |||
| 206 | static hash_varName_index *varName_ht = NULL; | ||
| 207 | |||
| 208 | ✗ | static void addVarToHashTable(const char *name, size_t index) | |
| 209 | { | ||
| 210 | hash_varName_index *s; | ||
| 211 | ✗ | HASH_FIND_STR(varName_ht, name, s); | |
| 212 | ✗ | assertStreamPrint(NULL, s, "Variable %s was not initialized in the DAG hash table.", name); | |
| 213 | ✗ | if (s->eqIndex != (size_t)(-1)) { | |
| 214 | ✗ | errorStreamPrint(OMC_LOG_STDOUT, 1, "Variable %s is solved in more than one equation.", name); | |
| 215 | ✗ | errorStreamPrint(OMC_LOG_STDOUT, 0, "originally solved in %zu, now in %zu", s->eqIndex, index); | |
| 216 | ✗ | messageClose(OMC_LOG_STDOUT); | |
| 217 | } else { | ||
| 218 | ✗ | s->eqIndex = index; | |
| 219 | } | ||
| 220 | ✗ | } | |
| 221 | |||
| 222 | ✗ | static void clearHashTable(void) | |
| 223 | { | ||
| 224 | hash_varName_index *s, *tmp; | ||
| 225 | |||
| 226 | ✗ | HASH_ITER(hh, varName_ht, s, tmp) { | |
| 227 | ✗ | HASH_DEL(varName_ht, s); | |
| 228 | ✗ | free(s); | |
| 229 | } | ||
| 230 | ✗ | varName_ht = NULL; | |
| 231 | ✗ | } | |
| 232 | |||
| 233 | ✗ | static void buildVarNameHashTable(MODEL_DATA *modelData) | |
| 234 | { | ||
| 235 | ✗ | for (int i = 0; i < modelData->nVariablesReal; ++i) { | |
| 236 | ✗ | const char *name = modelData->realVarsData[i].info.name; | |
| 237 | hash_varName_index *s; | ||
| 238 | ✗ | HASH_FIND_STR(varName_ht, name, s); | |
| 239 | ✗ | if (!s) { | |
| 240 | ✗ | s = (hash_varName_index *)malloc(sizeof *s); | |
| 241 | ✗ | s->id = name; | |
| 242 | ✗ | s->varIndex = modelData->realVarsData[i].info.id >= 1000 ? (size_t)(modelData->realVarsData[i].info.id - 1000) : (size_t)(-1); | |
| 243 | ✗ | s->eqIndex = (size_t)(-1); | |
| 244 | ✗ | HASH_ADD_KEYPTR(hh, varName_ht, s->id, strlen(s->id), s); | |
| 245 | } | ||
| 246 | } | ||
| 247 | ✗ | } | |
| 248 | |||
| 249 | /** | ||
| 250 | * @brief Get variable index from name | ||
| 251 | * | ||
| 252 | * @param name Name of variable | ||
| 253 | * @return size_t Index of variable if it exists else -1 | ||
| 254 | */ | ||
| 255 | ✗ | static size_t varIndexFromName(const char *name) | |
| 256 | { | ||
| 257 | hash_varName_index *s; | ||
| 258 | |||
| 259 | ✗ | HASH_FIND_STR(varName_ht, name, s); | |
| 260 | ✗ | return s ? s->varIndex : (size_t)(-1); | |
| 261 | } | ||
| 262 | |||
| 263 | /** | ||
| 264 | * @brief Set up DAG based on modelInfo | ||
| 265 | * | ||
| 266 | * Reads system info from modelData and sets all edges and var->eqn map in DAG. | ||
| 267 | * | ||
| 268 | * @param modelData Pointer to model data structure | ||
| 269 | * @param nEqns Number of equations in this DAG | ||
| 270 | * @param ixs Map from local eqIndices in this DAG to global eqIndices | ||
| 271 | */ | ||
| 272 | ✗ | void buildEvalDAG_ODE(MODEL_DATA *modelData, size_t nEqns, const size_t* ixs) | |
| 273 | { | ||
| 274 | size_t nEdges = 0; | ||
| 275 | ✗ | EVAL_DAG *dag = allocEvalDAG(modelData->nVariablesReal, nEqns); | |
| 276 | ✗ | modelData->dag = dag; | |
| 277 | |||
| 278 | ✗ | buildVarNameHashTable(modelData); | |
| 279 | |||
| 280 | ✗ | for (size_t i = 0; i < dag->nEqns; ++i) { | |
| 281 | ✗ | EQUATION_INFO eqInfo = modelInfoGetEquation(&modelData->modelDataXml, ixs[i]); | |
| 282 | |||
| 283 | /* see what variables it defines */ | ||
| 284 | ✗ | for (size_t j = 0; j < eqInfo.numVar; ++j) { | |
| 285 | ✗ | size_t varIndex = varIndexFromName(eqInfo.vars[j]); | |
| 286 | ✗ | if (varIndex != (size_t)(-1)) { | |
| 287 | /* store (varName -> eqIndex) in hash table */ | ||
| 288 | ✗ | addVarToHashTable(eqInfo.vars[j], i); | |
| 289 | |||
| 290 | /* set var -> eqn map */ | ||
| 291 | ✗ | dag->mapVarToEqNode[varIndex] = i; | |
| 292 | } | ||
| 293 | } | ||
| 294 | |||
| 295 | /* count how many edges will be in the DAG */ | ||
| 296 | // TODO remove duplicates if two vars are solved in the same eqn | ||
| 297 | ✗ | nEdges += eqInfo.numVarUsed; | |
| 298 | ✗ | dag->nEqDep[i] = eqInfo.numVarUsed; | |
| 299 | } | ||
| 300 | |||
| 301 | /* allocate space for all edges in one go */ | ||
| 302 | ✗ | dag->eqDep[0] = (size_t*) malloc(nEdges * sizeof(size_t)); | |
| 303 | ✗ | for (size_t i = 1; i < dag->nEqns; ++i) { | |
| 304 | ✗ | dag->eqDep[i] = dag->eqDep[i-1] + dag->nEqDep[i-1]; | |
| 305 | } | ||
| 306 | |||
| 307 | ✗ | for (size_t i = 0; i < dag->nEqns; ++i) { | |
| 308 | ✗ | EQUATION_INFO eqInfo = modelInfoGetEquation(&modelData->modelDataXml, ixs[i]); | |
| 309 | |||
| 310 | ✗ | for (size_t j = 0; j < eqInfo.numVarUsed; ++j) { | |
| 311 | /* look up variable in hash table */ | ||
| 312 | hash_varName_index *s; | ||
| 313 | ✗ | HASH_FIND_STR(varName_ht, eqInfo.varsUsed[j], s); | |
| 314 | /* if it exists, add the corresponding equation to the dependency */ | ||
| 315 | // TODO reduce size of eqDep if used variables are not solved in the | ||
| 316 | // system corresponding to this DAG | ||
| 317 | ✗ | dag->eqDep[i][j] = s ? s->eqIndex : (size_t)(-1); | |
| 318 | } | ||
| 319 | } | ||
| 320 | |||
| 321 | ✗ | clearHashTable(); | |
| 322 | ✗ | } | |
| 323 | |||
| 324 | /** | ||
| 325 | * @brief Set up DAG based on modelInfo | ||
| 326 | * | ||
| 327 | * Reads system info from modelData and sets all edges and var->eqn map in DAG. | ||
| 328 | * | ||
| 329 | * @param jacobian Pointer to the jacobian | ||
| 330 | * @param modelData Pointer to model data structure | ||
| 331 | * @param nEqns Number of equations in this DAG | ||
| 332 | * @param ixs Map from local eqIndices in this DAG to global eqIndices | ||
| 333 | */ | ||
| 334 | ✗ | void buildEvalDAG_Jac(JACOBIAN* jacobian, MODEL_DATA *modelData, size_t nEqns, const size_t* ixs) | |
| 335 | { | ||
| 336 | size_t nEdges = 0; | ||
| 337 | ✗ | EVAL_DAG *dag = allocEvalDAG(jacobian->sizeRows + jacobian->sizeTmpVars, nEqns); | |
| 338 | ✗ | jacobian->dag = dag; | |
| 339 | |||
| 340 | ✗ | buildVarNameHashTable(modelData); | |
| 341 | |||
| 342 | ✗ | for (size_t i = 0; i < dag->nEqns; ++i) { | |
| 343 | ✗ | EQUATION_INFO eqInfo = modelInfoGetEquation(&modelData->modelDataXml, ixs[i]); | |
| 344 | |||
| 345 | /* see what variables it defines */ | ||
| 346 | ✗ | for (size_t j = 0; j < eqInfo.numVar; ++j) { | |
| 347 | ✗ | size_t varIndex = varIndexFromName(eqInfo.vars[j]); | |
| 348 | ✗ | if (varIndex != (size_t)(-1)) { | |
| 349 | /* store (varName -> eqIndex) in hash table */ | ||
| 350 | ✗ | addVarToHashTable(eqInfo.vars[j], i); | |
| 351 | |||
| 352 | /* set var -> eqn map */ | ||
| 353 | ✗ | dag->mapVarToEqNode[varIndex] = i; | |
| 354 | } | ||
| 355 | } | ||
| 356 | |||
| 357 | /* count how many edges will be in the DAG */ | ||
| 358 | // TODO remove duplicates if two vars are solved in the same eqn | ||
| 359 | ✗ | nEdges += eqInfo.numVarUsed; | |
| 360 | ✗ | dag->nEqDep[i] = eqInfo.numVarUsed; | |
| 361 | } | ||
| 362 | |||
| 363 | /* allocate space for all edges in one go */ | ||
| 364 | ✗ | dag->eqDep[0] = (size_t*) malloc(nEdges * sizeof(size_t)); | |
| 365 | ✗ | for (size_t i = 1; i < dag->nEqns; ++i) { | |
| 366 | ✗ | dag->eqDep[i] = dag->eqDep[i-1] + dag->nEqDep[i-1]; | |
| 367 | } | ||
| 368 | |||
| 369 | ✗ | for (size_t i = 0; i < dag->nEqns; ++i) { | |
| 370 | ✗ | EQUATION_INFO eqInfo = modelInfoGetEquation(&modelData->modelDataXml, ixs[i]); | |
| 371 | |||
| 372 | ✗ | for (size_t j = 0; j < eqInfo.numVarUsed; ++j) { | |
| 373 | /* look up variable in hash table */ | ||
| 374 | hash_varName_index *s; | ||
| 375 | ✗ | HASH_FIND_STR(varName_ht, eqInfo.varsUsed[j], s); | |
| 376 | /* if it exists, add the corresponding equation to the dependency */ | ||
| 377 | // TODO reduce size of eqDep if used variables are not solved in the | ||
| 378 | // system corresponding to this DAG | ||
| 379 | ✗ | dag->eqDep[i][j] = s ? s->eqIndex : (size_t)(-1); | |
| 380 | } | ||
| 381 | } | ||
| 382 | |||
| 383 | ✗ | clearHashTable(); // done | |
| 384 | ✗ | } | |
| 385 |