diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c index ab526aea199..d15b1ca5824 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c @@ -1874,6 +1874,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn /* update time with performed stepSize */ gbData->time += gbData->lastStepSize; + gbData->timeDense = gbData->time; /* step is accepted and yOld needs to be updated */ memcpy(gbData->yOld, gbData->y, nStates * sizeof(double)); @@ -1924,25 +1925,24 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn // Current solution: Step back to the communication interval before the event and event detection // needs to be repeated listClear(solverInfo->eventLst); - gbData->lastStepSize = (eventTime - solverInfo->currentStepSize/2) - gbData->time; + gbData->lastStepSize = (eventTime - solverInfo->currentStepSize/2) - gbData->timeLeft; sData->timeValue = (eventTime - solverInfo->currentStepSize/2); gb_interpolation(gbData->interpolation, gbData->timeLeft, gbData->yLeft, gbData->kLeft, gbData->timeRight, gbData->yRight, gbData->kRight, sData->timeValue, sData->realVars, nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k); - memcpy(gbData->y, sData->realVars, gbData->nStates * sizeof(double)); - + memcpy(gbData->yOld, sData->realVars, gbData->nStates * sizeof(double)); gbData->timeRight = sData->timeValue; + gbData->time = gbData->timeRight; memcpy(gbData->yRight, sData->realVars, gbData->nStates * sizeof(double)); gbode_fODE(data, threadData, &(gbData->stats.nCallsODE)); memcpy(gbData->kRight, fODE, nStates * sizeof(double)); - } } infoStreamPrint(LOG_SOLVER, 0, "Accept step from %10g to %10g, error %10g interpolation error %10g, new stepsize %10g", - gbData->time - gbData->lastStepSize, gbData->time, err, gbData->err_int, gbData->stepSize); + gbData->timeLeft, gbData->timeRight, err, gbData->err_int, gbData->stepSize); /* emit step, if integratorSteps is selected */ if (solverInfo->integratorSteps) @@ -1980,11 +1980,23 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn // use chosen interpolation for emitting equidistant output (default hermite) if (solverInfo->currentStepSize > 0) { - gb_interpolation(gbData->interpolation, + if (gbData->timeDense > gbData->timeRight && (gbData->interpolation == GB_DENSE_OUTPUT || gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL)) + { + /* This case is needed, if an event has been detected during a large step (gbData->timeDense) of the integration + * and the integrator (gbData->timeRight) has been set back to the time just before the event. In this case the + * values in gbData->x and gbData->k are correct for the overall time intervall from gbData->timeLeft to gbData->timeDense */ + gb_interpolation(gbData->interpolation, + gbData->timeLeft, gbData->yLeft, gbData->kLeft, + gbData->timeDense, gbData->yRight, gbData->kRight, + sData->timeValue, sData->realVars, + nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k); + } else { + gb_interpolation(gbData->interpolation, gbData->timeLeft, gbData->yLeft, gbData->kLeft, gbData->timeRight, gbData->yRight, gbData->kRight, sData->timeValue, sData->realVars, nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k); + } } // log the emitted result if (ACTIVE_STREAM(LOG_GBODE)){ diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.h b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.h index 1a66a375e68..77f26b57286 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.h @@ -145,7 +145,7 @@ typedef struct DATA_GBODE{ double *stepSizeValues; /* ring buffer for step size control */ double err_slow, err_fast, err_int; /* error of the slow, fast states and a preiction of the interpolation error */ 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 time, timeLeft, timeRight, timeDense; /* actual time values and the time values of the current interpolation interval and for dense output */ double stepSize, lastStepSize, optStepSize; /* actual, last, and optimal step size of integration */ double maxStepSize; /* maximal step size of integration */ double initialStepSize; /* initial step size of integration */