Skip to content

Commit

Permalink
Utilize max and initial step size (#10334)
Browse files Browse the repository at this point in the history 
* Introduce allowed maximal step size
* Handling of flag initialStepSize added
  • Loading branch information
@bernhardbachmann
bernhardbachmann committed Mar 7, 2023
1 parent ca6dc41 commit 7a6d3ac
Show file tree
Hide file tree
Showing 3 changed files with 85 additions and 41 deletions.
85 changes: 45 additions & 40 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_ctrl.c
View file
Expand Up @@ -187,58 +187,63 @@ void getInitStepSize(DATA* data, threadData_t* threadData, DATA_GBODE* gbData)
gbData->time = sData->timeValue;
memcpy(gbData->yOld, sData->realVars, nStates*sizeof(double));
gbode_fODE(data, threadData, &(gbData->stats.nCallsODE));
memcpy(gbData->f, fODE, nStates*sizeof(double));

for (i=0; i<nStates; i++) {
sc = absTol + fabs(sDataOld->realVars[i])*relTol;
d0 += ((sDataOld->realVars[i] * sDataOld->realVars[i])/(sc*sc));
d1 += ((fODE[i] * fODE[i]) / (sc*sc));
}
d0 /= nStates;
d1 /= nStates;
if (gbData->initialStepSize < 0) {
memcpy(gbData->f, fODE, nStates*sizeof(double));
for (i=0; i<nStates; i++) {
sc = absTol + fabs(sDataOld->realVars[i])*relTol;
d0 += ((sDataOld->realVars[i] * sDataOld->realVars[i])/(sc*sc));
d1 += ((fODE[i] * fODE[i]) / (sc*sc));
}
d0 /= nStates;
d1 /= nStates;

d0 = sqrt(d0);
d1 = sqrt(d1);
d0 = sqrt(d0);
d1 = sqrt(d1);

/* calculate first guess of the initial step size */
if (d0 < 1e-5 || d1 < 1e-5) {
h0 = 1e-6;
} else {
h0 = 0.01 * d0/d1;
}
if (gbData->initialFailures>0)
h0 /= pow(10,gbData->initialFailures);
/* calculate first guess of the initial step size */
if (d0 < 1e-5 || d1 < 1e-5) {
h0 = 1e-6;
} else {
h0 = 0.01 * d0/d1;
}
if (gbData->initialFailures>0)
h0 /= pow(10,gbData->initialFailures);

for (i=0; i<nStates; i++) {
sData->realVars[i] = gbData->yOld[i] + fODE[i] * h0;
}
sData->timeValue += h0;
for (i=0; i<nStates; i++) {
sData->realVars[i] = gbData->yOld[i] + fODE[i] * h0;
}
sData->timeValue += h0;

gbode_fODE(data, threadData, &(gbData->stats.nCallsODE));
gbode_fODE(data, threadData, &(gbData->stats.nCallsODE));

for (i=0; i<nStates; i++) {
sc = absTol + fabs(gbData->yOld[i])*relTol;
d2 += ((fODE[i]-gbData->f[i])*(fODE[i]-gbData->f[i])/(sc*sc));
}
for (i=0; i<nStates; i++) {
sc = absTol + fabs(gbData->yOld[i])*relTol;
d2 += ((fODE[i]-gbData->f[i])*(fODE[i]-gbData->f[i])/(sc*sc));
}

d2 /= h0;
d2 = sqrt(d2);
d2 /= h0;
d2 = sqrt(d2);


d = fmax(d1,d2);
d = fmax(d1,d2);

if (d > 1e-15) {
h1 = sqrt(0.01/d);
} else {
h1 = fmax(1e-6, h0*1e-3);
}
if (d > 1e-15) {
h1 = sqrt(0.01/d);
} else {
h1 = fmax(1e-6, h0*1e-3);
}

gbData->stepSize = 0.5*fmin(100*h0,h1);
gbData->lastStepSize = 0.0;
gbData->stepSize = 0.5*fmin(100*h0,h1);
gbData->lastStepSize = 0.0;

sData->timeValue = gbData->time;
memcpy(sData->realVars, gbData->yOld, nStates*sizeof(double));
memcpy(fODE, gbData->f, nStates*sizeof(double));
sData->timeValue = gbData->time;
memcpy(sData->realVars, gbData->yOld, nStates*sizeof(double));
memcpy(fODE, gbData->f, nStates*sizeof(double));
} else {
gbData->stepSize = gbData->initialStepSize;
gbData->lastStepSize = 0.0;
}

infoStreamPrint(LOG_SOLVER, 0, "Initial step size = %e at time %g", gbData->stepSize, gbData->time);

Expand Down
39 changes: 38 additions & 1 deletion OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c
View file
Expand Up @@ -355,7 +355,36 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver

gbData->ctrl_method = getControllerMethod(FLAG_SR_CTRL);
gbData->stepSize_control = getControllFunc(gbData->ctrl_method);

/* define maximum step size gbode is allowed to go */
if (omc_flag[FLAG_MAX_STEP_SIZE])
{
gbData->maxStepSize = atof(omc_flagValue[FLAG_MAX_STEP_SIZE]);
if (gbData->maxStepSize < 0 || gbData->maxStepSize > DBL_MAX/2) {
throwStreamPrint(NULL, "maximum step size %g is not allowed", gbData->maxStepSize);
} else {
infoStreamPrint(LOG_SOLVER, 0, "maximum step size %g", gbData->maxStepSize);
}
}
else
{
gbData->maxStepSize = -1;
infoStreamPrint(LOG_SOLVER, 0, "maximum step size not set");
}
/* Initial step size */
if (omc_flag[FLAG_INITIAL_STEP_SIZE])
{
gbData->initialStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
if (gbData->initialStepSize < GB_MINIMAL_STEP_SIZE || gbData->initialStepSize > DBL_MAX/2) {
throwStreamPrint(NULL, "initial step size %g is not allowed, minimal step size is %g", gbData->initialStepSize, GB_MINIMAL_STEP_SIZE);
} else {
infoStreamPrint(LOG_SOLVER, 0, "initial step size %g", gbData->initialStepSize);
}
}
else
{
gbData->initialStepSize = -1; /* use default */
infoStreamPrint(LOG_SOLVER, 0, "initial step size not set");
}
gbData->isFirstStep = TRUE;

/* Allocate internal memory */
Expand Down Expand Up @@ -1290,6 +1319,8 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
// Call the step size control
gbData->lastStepSize = gbData->stepSize;
gbData->stepSize *= gbData->stepSize_control(gbData->errValues, gbData->stepSizeValues, gbData->tableau->error_order);
if (gbData->maxStepSize > 0 && gbData->maxStepSize < gbData->stepSize)
gbData->stepSize = gbData->maxStepSize;

// reject step, if error is too large
if ((err > 1 ) && gbData->ctrl_method != GB_CTRL_CNST) {
Expand Down Expand Up @@ -1340,6 +1371,8 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
if (gbData->err_int> err) {
gbData->errValues[0] = gbData->err_int;
gbData->stepSize = gbData->lastStepSize * gbData->stepSize_control(gbData->errValues, gbData->stepSizeValues, gbData->tableau->error_order);
if (gbData->maxStepSize > 0 && gbData->maxStepSize < gbData->stepSize)
gbData->stepSize = gbData->maxStepSize;
}
}
// reject step, if interpolaton error is too large
Expand Down Expand Up @@ -1722,6 +1755,8 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
// Call the step size control
gbData->lastStepSize = gbData->stepSize;
gbData->stepSize *= gbData->stepSize_control(gbData->errValues, gbData->stepSizeValues, gbData->tableau->error_order);
if (gbData->maxStepSize > 0 && gbData->maxStepSize < gbData->stepSize)
gbData->stepSize = gbData->maxStepSize;

// reject step, if error is too large
if ((err > 1) && gbData->ctrl_method != GB_CTRL_CNST) {
Expand Down Expand Up @@ -1759,6 +1794,8 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
if (gbData->err_int> err) {
gbData->errValues[0] = gbData->err_int;
gbData->stepSize = gbData->lastStepSize * gbData->stepSize_control(gbData->errValues, gbData->stepSizeValues, gbData->tableau->error_order);
if (gbData->maxStepSize > 0 && gbData->maxStepSize < gbData->stepSize)
gbData->stepSize = gbData->maxStepSize;
}
}
// reject step, if interpolaton error is too large
Expand Down
2 changes: 2 additions & 0 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.h
View file
Expand Up @@ -147,6 +147,8 @@ typedef struct DATA_GBODE{
double percentage, err_threshold; /* percentage of fast states and the corresponding error threshold */
double time, timeLeft, timeRight; /* actual time values and the time values of the current interpolation interval */
double stepSize, lastStepSize; /* actual and last step size of integration */
double maxStepSize; /* maximal step size of integration */
double initialStepSize; /* initial step size of integration */
int act_stage; /* Current stage of Runge-Kutta method. */
enum GB_CTRL_METHOD ctrl_method; /* Step size control algorithm */
int ringBufferSize; /* Buffer size for storing the error, stepSize and last values of states (yv) and their derivatives (kv) */
Expand Down

0 comments on commit 7a6d3ac

Please sign in to comment.