我正在使用 CVODES 来计算一个非常基本的方程的伴随灵敏度，该方程在 t=1 处具有解不连续性。在积分正向解决方案时，我在循环中使用CVodeF()in模式在每个间隔上积分，调用以在不连续处重新开始积分。这产生了正确的正向解决方案。然后我呼吁对伴随敏感性分析进行后向整合。CV_ONE_STEPCVodeReInit()CVodeB()
我的问题是关于前向集成的重新启动以及在后向集成期间如何处理它们。在我的程序开始时，我CVodeAdjInit()用调用Nd = 1，所以我相信我在每个集成步骤都保存了检查点。但是，在检查检查点并尝试CVodeGetAdjY()在 t=1 及其前后调用前向解决方案（带有 ）时，我没有在正确的时间看到正确的跳跃不连续性。
这些重启是否明确保存（在检查点方案或其他地方）？我是否应该调用其他 CVODES 函数来通知向后积分器这些重新启动？任何使用 CVODES 进行伴随灵敏度分析的一般指导（在存在正解不连续性的情况下）都将不胜感激。
我已经包含了说明这一点的示例代码。我们将 dy/dt = -0.05*y 从 t = 0 积分到 t = 2，在 t = 1 处施加 0.1 的脉冲。在这个例子中，我没有对伴随状态 lambda 做任何事情——主要目的是说明在后向积分过程中如何调用前向解 y。
#include <stdio.h>
#include <cvodes/cvodes.h>
#include <nvector/nvector_serial.h> /* access to serial N_Vector */
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
/* Number of equations in system (1 for this example) */
#define NEQ 1
/* Accessor macros */
#define Ith(v, i) NV_Ith_S(v,i) /* i-th vector component, i=0..NEQ-1 */
#define IJth(A,i,j) SM_ELEMENT_D(A,i,j) /* IJth numbers rows,cols 0..NEQ-1 */
/* Decay rate*/
#define R 0.05;
static int rhs(realtype t, N_Vector y, N_Vector ydot, void *user_data);
static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
static int rhs_adj(realtype t, N_Vector y, N_Vector lam, N_Vector lamdot, void *user_data);
static int Jac_adj(realtype t,
N_Vector y, N_Vector lam, N_Vector fyB,
SUNMatrix JB, void *user_data,
N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B);
int main() {
uint32_t maxord = 5; //BDF order
double abstol = 1e-8;
double reltol = 1e-8;
double abstol_adj = 1e-8;
double reltol_adj = 1e-8;
int adj_steps = 1; //number of integration steps between two consecutive checkpoints
int ncheck; //number of checkpoints
int indexB; //index for the backward problem
/* impulse at t=1*/
double impulse_time = 1.0;
double impulse_amount = 0.1;
/* integrate to t=2 */
double final_time = 2.0;
/* y = state */
N_Vector y = N_VNew_Serial(NEQ);
Ith(y, 0) = 1.0; //init condit.
/* lambda in adjoint equation. Needed for the backward integration, though will be ignored for this example */
N_Vector lam = N_VNew_Serial(NEQ);
Ith(lam, 0) = 0.0; //init condit.
/* initialize cvodes, set tolerances, etc. */
void *cvode_mem = CVodeCreate(CV_BDF);
CVodeSetMaxOrd(cvode_mem, maxord);
CVodeInit(cvode_mem, rhs, 0.0, y);
CVodeSStolerances(cvode_mem, reltol, abstol);
/* Create SUNMatrix and linear solver for the forward problem and set Jacobian fcn*/
SUNMatrix A = SUNDenseMatrix(NEQ, NEQ);
SUNLinearSolver LS = SUNLinSol_Dense(y, A);
CVodeSetLinearSolver(cvode_mem, LS, A);
CVodeSetJacFn(cvode_mem, Jac);
/* Inform the forward problem there will be an adjoint problem to solve as well */
CVodeAdjInit(cvode_mem, adj_steps, CV_HERMITE);
/* Initialization steps for adj. problem*/
CVodeCreateB(cvode_mem, CV_BDF, &indexB);
CVodeSetMaxOrdB(cvode_mem, indexB, maxord);
CVodeInitB(cvode_mem, indexB, rhs_adj, final_time, lam);
CVodeSStolerancesB(cvode_mem, indexB, reltol_adj, abstol_adj);
/* Create SUNMatrix and linear solver for the adj problem and attach */
SUNMatrix AB = SUNDenseMatrix(NEQ, NEQ);
SUNLinearSolver LSB = SUNLinSol_Dense(y, AB);
CVodeSetLinearSolverB(cvode_mem, indexB, LSB, AB);
CVodeSetJacFnB(cvode_mem, indexB, Jac_adj);
/* The forward integration */
realtype time; // updated by each integration step by CVodeF
double current_time = 0.0;
double goal_time = impulse_time; // we will first integrate to the time of the impulse
int impulse_applied = 0;
int init_needed = 0;
while (current_time < final_time) {
/* need to re-initialize cvodes after the jump */
if (init_needed) {
CVodeReInit(cvode_mem, current_time, y);
printf("Re-init after impulse\n");
init_needed = 0;
}
while (current_time < goal_time) {
/* main forward integration step */
CVodeF(cvode_mem, goal_time, y, &time, CV_ONE_STEP, &ncheck);
current_time = time;
printf("t = %10.8f, y = %10.8f | ncheck = %d\n", current_time, Ith(y, 0), ncheck);
}
/* apply impulse */
if (!impulse_applied && impulse_time <= current_time) {
current_time = impulse_time;
CVodeGetDky(cvode_mem, impulse_time, 0, y);
printf("****** Before impulse: t = %10.8f, y = %10.8f \n", current_time, Ith(y, 0));
Ith(y, 0) += impulse_amount;
printf("****** After impulse: t = %10.8f, y = %10.8f \n", current_time, Ith(y, 0));
init_needed = 1;
impulse_applied = 1;
goal_time = final_time;
}
}
/* Now, integrate backwards in time for the adjoint problem */
goal_time = 1.0;
printf("\nPerforming backward integration ...\n");
while (current_time > goal_time) {
/* main backward integration step */
CVodeB(cvode_mem, goal_time, CV_ONE_STEP);
/* need to call CVodeGetB to get the time of the last CVodeB step */
CVodeGetB(cvode_mem, indexB, &time, lam);
/* obtain interpolated forward solution at current time as well */
CVodeGetAdjY(cvode_mem, time, y);
printf("t = %10.8f, y = %10.8f \n", time, Ith(y, 0));
current_time = time;
}
printf("Around impulse: \n");
double times[5] = {1.002, 1.001, 1.0, 0.999, 0.998};
for (int i = 0; i < 5; i++){
CVodeGetAdjY(cvode_mem, times[i], y);
printf("t = %10.8f, y = %10.8f \n", times[i], Ith(y, 0));
}
/* cleanup */
N_VDestroy(y);
N_VDestroy(lam);
CVodeFree(&cvode_mem);
SUNLinSolFree(LS);
SUNLinSolFree(LSB);
SUNMatDestroy(A);
SUNMatDestroy(AB);
return 0;
}
static int rhs(realtype t, N_Vector y, N_Vector ydot, void *user_data) {
/* exp decay */
double r = R;
Ith(ydot, 0) = -r * Ith(y, 0);
return (0);
}
static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) {
/* Jacobian of rhs */
double r = R;
IJth(J, 0, 0) = -r;
return (0);
}
static int rhs_adj(realtype t, N_Vector y, N_Vector lam, N_Vector lamdot, void *user_data) {
/* RHS of adjoint problem
* Note: the adjoint problem is lam' = -(J^*)lam - g_y^*, where J is the Jacobian matrix of the main problem.
* For this example, we take g = 0 everywhere */
double r = R;
Ith(lamdot, 0) = r * Ith(lam, 0);
return (0);
}
static int Jac_adj(realtype t,
N_Vector y, N_Vector lam, N_Vector fyB,
SUNMatrix JB, void *user_data,
N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B) {
/* Jacobian of adjoint problem */
double r = R;
IJth(JB, 0, 0) = r;
return (0);
}
在前向积分期间，我们有输出（大约 t = 1）
...
t = 0.77864484, y = 0.96181582 | ncheck = 18
t = 0.94641606, y = 0.95378132 | ncheck = 19
t = 1.11418727, y = 0.94581393 | ncheck = 20
****** Before impulse: t = 1.00000000, y = 0.95122937
****** After impulse: t = 1.00000000, y = 1.05122937
Re-init after impulse
t = 1.00197548, y = 1.05112555 | ncheck = 21
t = 1.00395097, y = 1.05102173 | ncheck = 22
...
在后向阶段（这里，y 是通过 获得的CVodeGetAdjY），
...
t = 1.00934328, y = 1.05073841
t = 1.00395097, y = 1.05102173
t = 1.00197548, y = 1.05112555
t = 0.98222065, y = 0.95207536
在 t = 1.00197548 处召回的 y 值当时是正确的（这是在前向积分期间采取的脉冲之后的第一步），但是如果我随后在脉冲附近查询 y（再次使用CVodeGetAdjY）：
Around impulse:
t = 1.00200000, y = 0.95113425
t = 1.00100000, y = 0.95118181
t = 1.00000000, y = 0.95122937
t = 0.99900000, y = 0.95127693
t = 0.99800000, y = 0.95132450
冲动不明显。在 t = 1.0 时调用的 y 是前脉冲值。看起来好像，即使CVodeReInit()在前向积分期间在脉冲之后立即调用，在后向积分期间也没有“看到”后脉冲 y。此外，如果后向积分器需要检查点之间的任何步骤，则 1.00197548 和 t = 1.0 之间的插值 y 可能会关闭。
换句话说，我的问题是：在重新启动转发问题之后，有没有办法确保这样的重新启动被保存并可以从检查点数据访问？