Skip to content

Commit a3b5314

Browse files
AnHeuermannAnHeuermann
authored
Get NLS Jacobian from Homotopy solver (#10792)
* Getter function for homotopy solver Jacobian * Refactor omc vector norm * Use infinity norm in GBODE instead of re-doing it * Fixing include error
1 parent 28817ae commit a3b5314

File tree

10 files changed

+110
-67
lines changed

10 files changed

+110
-67
lines changed

‎OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@
5959
#include "model_help.h"
6060
#include "newtonIteration.h"
6161
#include "nonlinearSystem.h"
62+
#include "omc_math.h"
6263
#include "simulation/options.h"
6364
#include "simulation/results/simulation_result.h"
6465
#include "simulation/jacobian_util.h"
@@ -1276,12 +1277,11 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
12761277
}
12771278
}
12781279

1279-
// calculate corresponding values for error estimator and step size control (infinity norm)
1280-
for (i = 0, err=0; i < nStates; i++) {
1280+
// calculate corresponding values for error estimator and step size control
1281+
for (i = 0; i < nStates; i++) {
12811282
gbData->errtol[i] = Rtol * fmax(fabs(gbData->y[i]), fabs(gbData->yOld[i])) + Atol;
12821283
gbData->errest[i] = fabs(gbData->y[i] - gbData->yt[i]);
12831284
gbData->err[i] = gbData->tableau->fac * gbData->errest[i] / gbData->errtol[i];
1284-
err = fmax(err, gbData->err[i]);
12851285
}
12861286

12871287
if (ACTIVE_STREAM(LOG_GBODE_V))
@@ -1763,12 +1763,11 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
17631763
}
17641764
}
17651765

1766-
// calculate corresponding values for error estimator and step size control (infinity norm)
1767-
for (i = 0, err=0; i < nStates; i++) {
1766+
// calculate corresponding values for error estimator and step size control
1767+
for (i = 0; i < nStates; i++) {
17681768
gbData->errtol[i] = Rtol * fmax(fabs(gbData->y[i]), fabs(gbData->yOld[i])) + Atol;
17691769
gbData->errest[i] = fabs(gbData->y[i] - gbData->yt[i]);
17701770
gbData->err[i] = gbData->tableau->fac * gbData->errest[i] / gbData->errtol[i];
1771-
err = fmax(err, gbData->err[i]);
17721771
}
17731772

17741773
// Rotate buffer
@@ -1777,7 +1776,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
17771776
gbData->stepSizeValues[i] = gbData->stepSizeValues[i - 1];
17781777
}
17791778
// update new values
1780-
gbData->errValues[0] = err;
1779+
gbData->errValues[0] = _omc_gen_maximumVectorNorm(gbData->err, nStates);
17811780
gbData->stepSizeValues[0] = gbData->stepSize;
17821781

17831782
// Store performed step size for latter interpolation

‎OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -139,7 +139,8 @@ typedef struct DATA_GBODE{
139139
double *yv, *kv, *tv; /* Buffer storage of the last values of states (yv) and their derivatives (kv) */
140140
double *yr, *kr, *tr; /* Backup storage of the buffers yv, kv, tv for Richardson extrapolation */
141141
double *res_const; /* Constant parts of residual for non-linear system of implicit RK method. */
142-
double *errest, *errtol; /* absolute error and given error tolerance of each individual states */
142+
double *errest; /* Absolute error estimator of each individual state */
143+
double *errtol; /* Given error tolerance of each individual state */
143144
double *err; /* error of each individual state during integration err = errest/errtol*/
144145
double *errValues; /* ring buffer for step size control */
145146
double *stepSizeValues; /* ring buffer for step size control */

‎OMCompiler/SimulationRuntime/c/simulation/solver/newtonIteration.c‎

Lines changed: 35 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -40,9 +40,10 @@ extern "C" {
4040
#include <string.h> /* memcpy */
4141

4242
#include "simulation/simulation_info_json.h"
43+
#include "model_help.h"
44+
#include "omc_math.h"
4345
#include "util/omc_error.h"
4446
#include "util/varinfo.h"
45-
#include "model_help.h"
4647

4748
#include "nonlinearSystem.h"
4849
#include "newtonIteration.h"
@@ -456,24 +457,44 @@ void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x
456457
*scaledError_f = enorm_(&n,solverData->fvecScaled);
457458
}
458459

459-
/*! \fn calculatingErrors
460+
/**
461+
* @brief Compute residual scaling vector.
460462
*
461-
* function scales the residual vector using the jacobian (heuristic)
463+
* scalingVector[i] = 1 / ||Jac(i,:)||
464+
* Warn if Jacobian row is all zeros i.e. the Jacobian is singular.
465+
*
466+
* @param solverData Newton solver data.
467+
* @param scalingVector Residual scaling vector.
462468
*/
463-
void scaling_residual_vector(DATA_NEWTON* solverData)
464-
{
465-
int i,j,k;
466-
for(i=0, k=0; i<solverData->n; i++)
469+
void compute_scaling_vector(DATA_NEWTON* solverData, double* scalingVector) {
470+
int i;
471+
int jac_row_start;
472+
473+
for(i=0; i<solverData->n; i++)
467474
{
468-
solverData->resScaling[i] = 0.0;
469-
for(j=0; j<solverData->n; j++, ++k)
470-
{
471-
solverData->resScaling[i] = fmax(fabs(solverData->fjac[k]), solverData->resScaling[i]);
472-
}
473-
if(solverData->resScaling[i] <= 0.0){
475+
jac_row_start = i*solverData->n;
476+
scalingVector[i] = _omc_gen_maximumVectorNorm(&(solverData->fjac[jac_row_start]), solverData->n);
477+
if(scalingVector[i] <= 0.0) {
474478
warningStreamPrint(LOG_NLS_V, 1, "Jacobian matrix is singular.");
475-
solverData->resScaling[i] = 1e-16;
479+
scalingVector[i] = 1e-16;
476480
}
481+
}
482+
}
483+
484+
/**
485+
* @brief Scale residual vector.
486+
*
487+
* Save result in solverData->fvecScaled.
488+
*
489+
* @param solverData Newton solver data.
490+
*/
491+
void scaling_residual_vector(DATA_NEWTON* solverData)
492+
{
493+
int i;
494+
495+
compute_scaling_vector(solverData, solverData->resScaling);
496+
for(i=0; i<solverData->n; i++)
497+
{
477498
solverData->fvecScaled[i] = solverData->fvec[i] / solverData->resScaling[i];
478499
}
479500
}

‎OMCompiler/SimulationRuntime/c/simulation/solver/newtonIteration.h‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -47,13 +47,13 @@ extern "C" {
4747
*/
4848
typedef struct DATA_NEWTON
4949
{
50-
int initialized; /**< 1 = initialized, else = 0*/
51-
double* resScaling;
52-
double* fvecScaled;
50+
int initialized; /** 1 = initialized, else = 0 */
51+
double* resScaling; /** Residual scaling vector */
52+
double* fvecScaled; /** scaled f vector */
5353

5454
int newtonStrategy;
5555

56-
int n; /**< size of equation */
56+
int n; /** size of equation system */
5757
double* x;
5858
double* fvec;
5959
double xtol;
@@ -62,13 +62,13 @@ typedef struct DATA_NEWTON
6262
int maxfev;
6363
int info;
6464
double epsfcn;
65-
double* fjac;
65+
double* fjac; /** Jacobian matrix in row-major format */
6666
double* rwork;
6767
int* iwork;
6868
int calculate_jacobian;
6969
int factorization;
70-
int numberOfIterations; /**< over the whole simulation time */
71-
int numberOfFunctionEvaluations; /**< over the whole simulation time */
70+
int numberOfIterations; /** over the whole simulation time */
71+
int numberOfFunctionEvaluations; /** over the whole simulation time */
7272

7373
/* damped newton */
7474
double* x_new;
@@ -81,8 +81,8 @@ typedef struct DATA_NEWTON
8181
rtclock_t timeClock;
8282

8383
/* Debug information */
84-
double time; /**< Simulation time */
85-
modelica_boolean initial; /**< True if in initialization */
84+
double time; /** Simulation time */
85+
modelica_boolean initial; /** True if in initialization */
8686

8787
NLS_USERDATA* userData;
8888
} DATA_NEWTON;

‎OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -114,7 +114,7 @@ typedef struct DATA_HOMOTOPY
114114
double* x1;
115115
double* finit;
116116
double* fx0;
117-
double* fJac;
117+
double* fJac; /* n times n Jacobian matrix with additional scaling row at the end */
118118
double* fJacx0;
119119

120120
/* debug arrays */
@@ -144,7 +144,7 @@ typedef struct DATA_HOMOTOPY
144144

145145
int (*f) (struct DATA_HOMOTOPY*, double*, double*);
146146
int (*f_con) (struct DATA_HOMOTOPY*, double*, double*);
147-
int (*fJac_f) (struct DATA_HOMOTOPY*, double*, double*);
147+
int (*fJac_f) (struct DATA_HOMOTOPY* solverData, double* x, double* fJac);
148148
int (*h_function)(struct DATA_HOMOTOPY*, double*, double*);
149149
int (*hJac_dh) (struct DATA_HOMOTOPY*, double*, double*);
150150

@@ -2532,4 +2532,15 @@ NLS_SOLVER_STATUS solveHomotopy(DATA *data, threadData_t *threadData, NONLINEAR_
25322532
return success;
25332533
}
25342534

2535+
/**
2536+
* @brief Return pointer to Jacobian.
2537+
*
2538+
* @param nlsData Non-linear system data.
2539+
* @return double* Jacobian in row-major format.
2540+
*/
2541+
double* getHomotopyJacobian(NONLINEAR_SYSTEM_DATA* nlsData) {
2542+
DATA_HOMOTOPY* homotopyData = (DATA_HOMOTOPY*)(nlsData->solverData);
2543+
return homotopyData->fJac;
2544+
}
2545+
25352546
#endif

‎OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@ typedef struct DATA_HOMOTOPY DATA_HOMOTOPY;
4747
DATA_HOMOTOPY* allocateHomotopyData(size_t size, NLS_USERDATA* userData);
4848
void freeHomotopyData(DATA_HOMOTOPY* homotopyData);
4949
NLS_SOLVER_STATUS solveHomotopy(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nlsData);
50+
double* getHomotopyJacobian(NONLINEAR_SYSTEM_DATA* nlsData);
5051

5152
#ifdef __cplusplus
5253
};

‎OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ typedef struct DATA_HYBRD
6969
integer info;
7070
integer nfev;
7171
integer njev;
72-
double* fjac;
72+
double* fjac; /* Jacobian matrix in row-major order */
7373
double* fjacobian;
7474
integer ldfjac;
7575
double* r__;

‎OMCompiler/SimulationRuntime/c/simulation/solver/omc_math.c‎

Lines changed: 40 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -783,63 +783,73 @@ void _omc_printMatrix(_omc_matrix* mat, const char* name, const enum LOG_STREAM
783783
}
784784
}
785785

786-
/*! \fn _omc_scalar _omc_euclideanVectorNorm(_omc_vector* vec)
786+
/**
787+
* @brief Euclidean vector norm.
787788
*
788-
* calculates the euclidean vector norm
789+
* ||x||_2 := sqrt(x_1^2 + ... + x_n^2)
789790
*
790-
* \param [in] [vec] !TODO: DESCRIBE ME!
791+
* @param vec OMC vector.
792+
* @return _omc_scalar Vector norm.
791793
*/
792794
_omc_scalar _omc_euclideanVectorNorm(const _omc_vector* vec)
793795
{
794-
_omc_size i;
795-
_omc_scalar result = 0;
796-
assertStreamPrint(NULL, vec->size > 0, "Vector size is greater than zero");
797-
assertStreamPrint(NULL, NULL != vec->data, "Vector data is NULL pointer");
798-
for (i = 0; i < vec->size; ++i) {
799-
result += pow(fabs(vec->data[i]),2.0);
800-
}
801-
802-
return sqrt(result);
796+
return _omc_gen_euclideanVectorNorm(vec->data, vec->size);
803797
}
804798

805-
/*! \fn _omc_scalar _omc_gen_euclideanVectorNorm(_omc_vector* vec)
799+
/**
800+
* @brief Euclidean vector norm.
806801
*
807-
* calculates the euclidean vector norm
802+
* ||x||_2 := sqrt(x_1^2 + ... + x_n^2)
808803
*
809-
* \param [in] [vec] !TODO: DESCRIBE ME!
804+
* @param vec_data Array with vector.
805+
* @param vec_size Length of array vec_data.
806+
* @return _omc_scalar Vector norm.
810807
*/
811808
_omc_scalar _omc_gen_euclideanVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size)
812809
{
813810
_omc_size i;
814811
_omc_scalar result = 0;
815812
assertStreamPrint(NULL, vec_size > 0, "Vector size is greater than zero");
816813
assertStreamPrint(NULL, NULL != vec_data, "Vector data is NULL pointer");
817-
for (i = 0; i < vec_size; ++i) {
818-
result += pow(fabs(vec_data[i]),2.0);
814+
for (i = 0; i < vec_size; i++) {
815+
result += vec_data[i]*vec_data[i];
819816
}
820817

821818
return sqrt(result);
822819
}
823820

824-
/*! \fn _omc_scalar _omc_maximumVectorNorm(_omc_vector* vec)
821+
/**
822+
* @brief Maximum / Infinity vector norm.
825823
*
826-
* calculates the maximum vector norm
824+
* ||x||_max := max(|x_1|, ..., |x_n|)
825+
*
826+
* @param vec OMC vector.
827+
* @return _omc_scalar Vector norm.
827828
*/
828829
_omc_scalar _omc_maximumVectorNorm(const _omc_vector* vec)
829830
{
831+
return _omc_gen_maximumVectorNorm(vec->data, vec->size);
832+
}
833+
834+
/**
835+
* @brief Maximum / Infinity vector norm.
836+
*
837+
* ||x||_max := max(|x_1|, ..., |x_n|)
838+
*
839+
* @param vec_data Array with vector.
840+
* @param vec_size Length of array vec_data.
841+
* @return _omc_scalar Vector norm.
842+
*/
843+
_omc_scalar _omc_gen_maximumVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size) {
830844
_omc_size i;
831-
_omc_scalar result, tmp;
832-
assertStreamPrint(NULL, vec->size > 0, "Vector size is greater the zero");
833-
assertStreamPrint(NULL, NULL != vec->data, "Vector data is NULL pointer");
834-
result = fabs(vec->data[0]);
835-
for (i = 1; i < vec->size; ++i)
845+
_omc_scalar norm;
846+
assertStreamPrint(NULL, vec_size > 0, "Vector size is greater the zero");
847+
assertStreamPrint(NULL, NULL != vec_data, "Vector data is NULL pointer");
848+
norm = fabs(vec_data[0]);
849+
for (i = 1; i < vec_size; i++)
836850
{
837-
tmp = fabs(vec->data[i]);
838-
if (result < tmp)
839-
{
840-
result = tmp;
841-
}
851+
norm = fmax(fabs(vec_data[i]), norm);
842852
}
843853

844-
return result;
854+
return norm;
845855
}

‎OMCompiler/SimulationRuntime/c/simulation/solver/omc_math.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -121,5 +121,6 @@ void _omc_printMatrix(_omc_matrix* mat, const char* name, const enum LOG_STREAM
121121
_omc_scalar _omc_euclideanVectorNorm(const _omc_vector* vec);
122122
_omc_scalar _omc_gen_euclideanVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size);
123123
_omc_scalar _omc_maximumVectorNorm(const _omc_vector* vec);
124+
_omc_scalar _omc_gen_maximumVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size);
124125

125126
#endif

‎OMCompiler/SimulationRuntime/c/simulation_data.h‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -781,10 +781,9 @@ typedef struct SIMULATION_INFO
781781

782782
ANALYTIC_JACOBIAN* analyticJacobians; // TODO Only store information for Jacobian used by integrator here
783783

784-
NONLINEAR_SYSTEM_DATA* nonlinearSystemData;
784+
NONLINEAR_SYSTEM_DATA* nonlinearSystemData; /* Array of non-linear systems */
785785

786-
LINEAR_SYSTEM_DATA* linearSystemData;
787-
int currentLinearSystemIndex;
786+
LINEAR_SYSTEM_DATA* linearSystemData; /* Array of linear systems */
788787

789788
MIXED_SYSTEM_DATA* mixedSystemData;
790789

0 commit comments

Comments
 (0)