Lines Matching refs:ts
45 static PetscErrorCode TSThetaGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot) in TSThetaGetX0AndXdot() argument
47 TS_Theta *th = (TS_Theta *)ts->data; in TSThetaGetX0AndXdot()
51 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSTheta_X0", X0)); in TSThetaGetX0AndXdot()
52 else *X0 = ts->vec_sol; in TSThetaGetX0AndXdot()
55 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSTheta_Xdot", Xdot)); in TSThetaGetX0AndXdot()
61 static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot) in TSThetaRestoreX0AndXdot() argument
65 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSTheta_X0", X0)); in TSThetaRestoreX0AndXdot()
68 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSTheta_Xdot", Xdot)); in TSThetaRestoreX0AndXdot()
81 TS ts = (TS)ctx; in DMRestrictHook_TSTheta() local
85 PetscCall(TSThetaGetX0AndXdot(ts, fine, &X0, &Xdot)); in DMRestrictHook_TSTheta()
86 PetscCall(TSThetaGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c)); in DMRestrictHook_TSTheta()
91 PetscCall(TSThetaRestoreX0AndXdot(ts, fine, &X0, &Xdot)); in DMRestrictHook_TSTheta()
92 PetscCall(TSThetaRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c)); in DMRestrictHook_TSTheta()
104 TS ts = (TS)ctx; in DMSubDomainRestrictHook_TSTheta() local
108 PetscCall(TSThetaGetX0AndXdot(ts, dm, &X0, &Xdot)); in DMSubDomainRestrictHook_TSTheta()
109 PetscCall(TSThetaGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub)); in DMSubDomainRestrictHook_TSTheta()
117 PetscCall(TSThetaRestoreX0AndXdot(ts, dm, &X0, &Xdot)); in DMSubDomainRestrictHook_TSTheta()
118 PetscCall(TSThetaRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub)); in DMSubDomainRestrictHook_TSTheta()
122 static PetscErrorCode TSThetaEvaluateCostIntegral(TS ts) in TSThetaEvaluateCostIntegral() argument
124 TS_Theta *th = (TS_Theta *)ts->data; in TSThetaEvaluateCostIntegral()
125 TS quadts = ts->quadraturets; in TSThetaEvaluateCostIntegral()
131 PetscCall(TSComputeRHSFunction(quadts, th->ptime0, th->X0, ts->vec_costintegrand)); in TSThetaEvaluateCostIntegral()
132 … PetscCall(VecAXPY(quadts->vec_sol, th->time_step0 * (1.0 - th->Theta), ts->vec_costintegrand)); in TSThetaEvaluateCostIntegral()
134 PetscCall(TSComputeRHSFunction(quadts, ts->ptime, ts->vec_sol, ts->vec_costintegrand)); in TSThetaEvaluateCostIntegral()
135 PetscCall(VecAXPY(quadts->vec_sol, th->time_step0 * th->Theta, ts->vec_costintegrand)); in TSThetaEvaluateCostIntegral()
137 PetscCall(TSComputeRHSFunction(quadts, th->stage_time, th->X, ts->vec_costintegrand)); in TSThetaEvaluateCostIntegral()
138 PetscCall(VecAXPY(quadts->vec_sol, th->time_step0, ts->vec_costintegrand)); in TSThetaEvaluateCostIntegral()
143 static PetscErrorCode TSForwardCostIntegral_Theta(TS ts) in TSForwardCostIntegral_Theta() argument
145 TS_Theta *th = (TS_Theta *)ts->data; in TSForwardCostIntegral_Theta()
146 TS quadts = ts->quadraturets; in TSForwardCostIntegral_Theta()
151 PetscCall(TSThetaEvaluateCostIntegral(ts)); in TSForwardCostIntegral_Theta()
155 static PetscErrorCode TSAdjointCostIntegral_Theta(TS ts) in TSAdjointCostIntegral_Theta() argument
157 TS_Theta *th = (TS_Theta *)ts->data; in TSAdjointCostIntegral_Theta()
161 th->ptime0 = ts->ptime + ts->time_step; in TSAdjointCostIntegral_Theta()
162 th->time_step0 = -ts->time_step; in TSAdjointCostIntegral_Theta()
163 PetscCall(TSThetaEvaluateCostIntegral(ts)); in TSAdjointCostIntegral_Theta()
167 static PetscErrorCode TSTheta_SNESSolve(TS ts, Vec b, Vec x) in TSTheta_SNESSolve() argument
172 PetscCall(SNESSolve(ts->snes, b, x)); in TSTheta_SNESSolve()
173 PetscCall(SNESGetIterationNumber(ts->snes, &nits)); in TSTheta_SNESSolve()
174 PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits)); in TSTheta_SNESSolve()
175 ts->snes_its += nits; in TSTheta_SNESSolve()
176 ts->ksp_its += lits; in TSTheta_SNESSolve()
180 static PetscErrorCode TSResizeRegister_Theta(TS ts, PetscBool reg) in TSResizeRegister_Theta() argument
182 TS_Theta *th = (TS_Theta *)ts->data; in TSResizeRegister_Theta()
186 PetscCall(TSResizeRegisterVec(ts, "ts:theta:sol_prev", th->vec_sol_prev)); in TSResizeRegister_Theta()
187 PetscCall(TSResizeRegisterVec(ts, "ts:theta:X0", th->X0)); in TSResizeRegister_Theta()
189 PetscCall(TSResizeRetrieveVec(ts, "ts:theta:sol_prev", &th->vec_sol_prev)); in TSResizeRegister_Theta()
191 PetscCall(TSResizeRetrieveVec(ts, "ts:theta:X0", &th->X0)); in TSResizeRegister_Theta()
197 static PetscErrorCode TSStep_Theta(TS ts) in TSStep_Theta() argument
199 TS_Theta *th = (TS_Theta *)ts->data; in TSStep_Theta()
202 PetscReal next_time_step = ts->time_step; in TSStep_Theta()
205 if (!ts->steprollback) { in TSStep_Theta()
207 PetscCall(VecCopy(ts->vec_sol, th->X0)); in TSStep_Theta()
211 while (!ts->reason && th->status != TS_STEP_COMPLETE) { in TSStep_Theta()
212 th->shift = 1 / (th->Theta * ts->time_step); in TSStep_Theta()
213 th->stage_time = ts->ptime + (th->endpoint ? (PetscReal)1 : th->Theta) * ts->time_step; in TSStep_Theta()
215 if (th->extrapolate && !ts->steprestart) PetscCall(VecAXPY(th->X, 1 / th->shift, th->Xdot)); in TSStep_Theta()
217 if (!th->affine) PetscCall(VecDuplicate(ts->vec_sol, &th->affine)); in TSStep_Theta()
219 PetscCall(TSComputeIFunction(ts, ts->ptime, th->X0, th->Xdot, th->affine, PETSC_FALSE)); in TSStep_Theta()
222 PetscCall(TSPreStage(ts, th->stage_time)); in TSStep_Theta()
223 PetscCall(TSTheta_SNESSolve(ts, th->endpoint ? th->affine : NULL, th->X)); in TSStep_Theta()
224 PetscCall(TSPostStage(ts, th->stage_time, 0, &th->X)); in TSStep_Theta()
225 PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->X, &stageok)); in TSStep_Theta()
229 PetscCall(VecCopy(th->X, ts->vec_sol)); in TSStep_Theta()
232 …if (th->Theta == 1.0) PetscCall(VecCopy(th->X, ts->vec_sol)); /* BEULER, stage alread… in TSStep_Theta()
234 PetscCall(VecAXPY(ts->vec_sol, ts->time_step, th->Xdot)); in TSStep_Theta()
235 … PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, ts->vec_sol, &stageok)); in TSStep_Theta()
237 PetscCall(VecCopy(th->X0, ts->vec_sol)); in TSStep_Theta()
244 PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); in TSStep_Theta()
247 PetscCall(VecCopy(th->X0, ts->vec_sol)); in TSStep_Theta()
248 ts->time_step = next_time_step; in TSStep_Theta()
252 …if (ts->forward_solve || ts->costintegralfwd) { /* Save the info for the later use in cost integra… in TSStep_Theta()
253 th->ptime0 = ts->ptime; in TSStep_Theta()
254 th->time_step0 = ts->time_step; in TSStep_Theta()
256 ts->ptime += ts->time_step; in TSStep_Theta()
257 ts->time_step = next_time_step; in TSStep_Theta()
261 ts->reject++; in TSStep_Theta()
263 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { in TSStep_Theta()
264 ts->reason = TS_DIVERGED_STEP_REJECTED; in TSStep_Theta()
265 …tscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than curr… in TSStep_Theta()
271 static PetscErrorCode TSAdjointStepBEuler_Private(TS ts) in TSAdjointStepBEuler_Private() argument
273 TS_Theta *th = (TS_Theta *)ts->data; in TSAdjointStepBEuler_Private()
274 TS quadts = ts->quadraturets; in TSAdjointStepBEuler_Private()
286 PetscCall(TSGetEquationType(ts, &eqtype)); in TSAdjointStepBEuler_Private()
289 VecsDeltaLam = ts->vecs_sensi; in TSAdjointStepBEuler_Private()
290 VecsDeltaLam2 = ts->vecs_sensi2; in TSAdjointStepBEuler_Private()
293 PetscCall(SNESGetKSP(ts->snes, &ksp)); in TSAdjointStepBEuler_Private()
294 PetscCall(TSGetIJacobian(ts, &J, &Jpre, NULL, NULL)); in TSAdjointStepBEuler_Private()
300 th->stage_time = ts->ptime; in TSAdjointStepBEuler_Private()
301 adjoint_time_step = -ts->time_step; /* always positive since time_step is negative */ in TSAdjointStepBEuler_Private()
304 if (quadts) PetscCall(TSComputeRHSJacobian(quadts, th->stage_time, ts->vec_sol, quadJ, NULL)); in TSAdjointStepBEuler_Private()
306 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStepBEuler_Private()
307 PetscCall(VecCopy(ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); in TSAdjointStepBEuler_Private()
311 PetscCall(VecPlaceArray(ts->vec_drdu_col, xarr)); in TSAdjointStepBEuler_Private()
312 PetscCall(VecAXPY(VecsSensiTemp[nadj], 1., ts->vec_drdu_col)); in TSAdjointStepBEuler_Private()
313 PetscCall(VecResetArray(ts->vec_drdu_col)); in TSAdjointStepBEuler_Private()
320 PetscCall(TSComputeSNESJacobian(ts, ts->vec_sol, J, Jpre)); in TSAdjointStepBEuler_Private()
324 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStepBEuler_Private()
329 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; in TSAdjointStepBEuler_Private()
330 …etscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve f… in TSAdjointStepBEuler_Private()
334 if (ts->vecs_sensi2) { /* U_{n+1} */ in TSAdjointStepBEuler_Private()
336 PetscCall(MatDenseGetColumn(ts->mat_sensip, 0, &xarr)); in TSAdjointStepBEuler_Private()
337 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); in TSAdjointStepBEuler_Private()
339 …scCall(TSComputeIHessianProductFunctionUU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_s… in TSAdjointStepBEuler_Private()
341 …PetscCall(TSComputeIHessianProductFunctionUP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->ve… in TSAdjointStepBEuler_Private()
342 for (nadj = 0; nadj < ts->numcost; nadj++) { /* compute the residual */ in TSAdjointStepBEuler_Private()
343 PetscCall(VecCopy(ts->vecs_sensi2[nadj], VecsSensi2Temp[nadj])); in TSAdjointStepBEuler_Private()
345 PetscCall(VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fuu[nadj])); in TSAdjointStepBEuler_Private()
346 if (ts->vecs_fup) PetscCall(VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fup[nadj])); in TSAdjointStepBEuler_Private()
349 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStepBEuler_Private()
354 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; in TSAdjointStepBEuler_Private()
355 …etscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve f… in TSAdjointStepBEuler_Private()
363 PetscCall(TSComputeSNESJacobian(ts, ts->vec_sol, J, Jpre)); in TSAdjointStepBEuler_Private()
365 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStepBEuler_Private()
370 PetscCall(VecCopy(VecsSensiTemp[nadj], ts->vecs_sensi[nadj])); in TSAdjointStepBEuler_Private()
371 if (ts->vecs_sensi2) { in TSAdjointStepBEuler_Private()
374 PetscCall(VecCopy(VecsSensi2Temp[nadj], ts->vecs_sensi2[nadj])); in TSAdjointStepBEuler_Private()
378 if (ts->vecs_sensip) { in TSAdjointStepBEuler_Private()
379 …PetscCall(VecAXPBYPCZ(th->Xdot, -1. / adjoint_time_step, 1.0 / adjoint_time_step, 0, th->X0, ts->v… in TSAdjointStepBEuler_Private()
380 …PetscCall(TSComputeIJacobianP(ts, th->stage_time, ts->vec_sol, th->Xdot, 1. / adjoint_time_step, t… in TSAdjointStepBEuler_Private()
381 if (quadts) PetscCall(TSComputeRHSJacobianP(quadts, th->stage_time, ts->vec_sol, quadJp)); in TSAdjointStepBEuler_Private()
382 if (ts->vecs_sensi2p) { in TSAdjointStepBEuler_Private()
384 …scCall(TSComputeIHessianProductFunctionPU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_s… in TSAdjointStepBEuler_Private()
386 …PetscCall(TSComputeIHessianProductFunctionPP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->ve… in TSAdjointStepBEuler_Private()
389 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStepBEuler_Private()
390 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj])); in TSAdjointStepBEuler_Private()
391 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsDeltaMu[nadj])); in TSAdjointStepBEuler_Private()
394 PetscCall(VecPlaceArray(ts->vec_drdp_col, xarr)); in TSAdjointStepBEuler_Private()
395 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, ts->vec_drdp_col)); in TSAdjointStepBEuler_Private()
396 PetscCall(VecResetArray(ts->vec_drdp_col)); in TSAdjointStepBEuler_Private()
399 if (ts->vecs_sensi2p) { in TSAdjointStepBEuler_Private()
400 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam2[nadj], VecsDeltaMu2[nadj])); in TSAdjointStepBEuler_Private()
401 PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step, VecsDeltaMu2[nadj])); in TSAdjointStepBEuler_Private()
402 …if (ts->vecs_fpu) PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step, ts->vecs_fpu[nadj]… in TSAdjointStepBEuler_Private()
403 …if (ts->vecs_fpp) PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step, ts->vecs_fpp[nadj]… in TSAdjointStepBEuler_Private()
408 if (ts->vecs_sensi2) { in TSAdjointStepBEuler_Private()
409 PetscCall(VecResetArray(ts->vec_sensip_col)); in TSAdjointStepBEuler_Private()
410 PetscCall(MatDenseRestoreColumn(ts->mat_sensip, &xarr)); in TSAdjointStepBEuler_Private()
416 static PetscErrorCode TSAdjointStep_Theta(TS ts) in TSAdjointStep_Theta() argument
418 TS_Theta *th = (TS_Theta *)ts->data; in TSAdjointStep_Theta()
419 TS quadts = ts->quadraturets; in TSAdjointStep_Theta()
431 PetscCall(TSAdjointStepBEuler_Private(ts)); in TSAdjointStep_Theta()
435 PetscCall(SNESGetKSP(ts->snes, &ksp)); in TSAdjointStep_Theta()
436 PetscCall(TSGetIJacobian(ts, &J, &Jpre, NULL, NULL)); in TSAdjointStep_Theta()
442 th->stage_time = th->endpoint ? ts->ptime : (ts->ptime + (1. - th->Theta) * ts->time_step); in TSAdjointStep_Theta()
443 adjoint_ptime = ts->ptime + ts->time_step; in TSAdjointStep_Theta()
444 adjoint_time_step = -ts->time_step; /* always positive since time_step is negative */ in TSAdjointStep_Theta()
448 …ecAXPBYPCZ(th->X0, 1.0 / (1.0 - th->Theta), th->Theta / (th->Theta - 1.0), 0, th->X, ts->vec_sol)); in TSAdjointStep_Theta()
455 PetscCall(TSComputeRHSJacobian(quadts, th->stage_time, ts->vec_sol, quadJ, NULL)); in TSAdjointStep_Theta()
461 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
462 PetscCall(VecCopy(ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); in TSAdjointStep_Theta()
466 PetscCall(VecPlaceArray(ts->vec_drdu_col, xarr)); in TSAdjointStep_Theta()
467 PetscCall(VecAXPY(VecsSensiTemp[nadj], 1., ts->vec_drdu_col)); in TSAdjointStep_Theta()
468 PetscCall(VecResetArray(ts->vec_drdu_col)); in TSAdjointStep_Theta()
476 PetscCall(TSComputeSNESJacobian(ts, ts->vec_sol, J, Jpre)); in TSAdjointStep_Theta()
478 PetscCall(TSComputeSNESJacobian(ts, th->X, J, Jpre)); in TSAdjointStep_Theta()
483 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
488 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; in TSAdjointStep_Theta()
489 …etscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve f… in TSAdjointStep_Theta()
494 if (ts->vecs_sensi2) { /* U_{n+1} */ in TSAdjointStep_Theta()
495 …PetscCheck(th->endpoint, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Operation not implement… in TSAdjointStep_Theta()
497 PetscCall(MatDenseGetColumn(ts->mat_sensip, 0, &xarr)); in TSAdjointStep_Theta()
498 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); in TSAdjointStep_Theta()
500 …scCall(TSComputeIHessianProductFunctionUU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_s… in TSAdjointStep_Theta()
501 PetscCall(VecResetArray(ts->vec_sensip_col)); in TSAdjointStep_Theta()
502 PetscCall(MatDenseRestoreColumn(ts->mat_sensip, &xarr)); in TSAdjointStep_Theta()
504 …PetscCall(TSComputeIHessianProductFunctionUP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->ve… in TSAdjointStep_Theta()
505 for (nadj = 0; nadj < ts->numcost; nadj++) { /* compute the residual */ in TSAdjointStep_Theta()
506 PetscCall(VecCopy(ts->vecs_sensi2[nadj], VecsSensi2Temp[nadj])); in TSAdjointStep_Theta()
508 PetscCall(VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fuu[nadj])); in TSAdjointStep_Theta()
509 if (ts->vecs_fup) PetscCall(VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fup[nadj])); in TSAdjointStep_Theta()
512 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
517 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; in TSAdjointStep_Theta()
518 …etscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve f… in TSAdjointStep_Theta()
527 PetscCall(TSComputeSNESJacobian(ts, th->X0, J, Jpre)); in TSAdjointStep_Theta()
531 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
532 PetscCall(MatMultTranspose(J, VecsDeltaLam[nadj], ts->vecs_sensi[nadj])); in TSAdjointStep_Theta()
535 PetscCall(VecPlaceArray(ts->vec_drdu_col, xarr)); in TSAdjointStep_Theta()
536 PetscCall(VecAXPY(ts->vecs_sensi[nadj], -1., ts->vec_drdu_col)); in TSAdjointStep_Theta()
537 PetscCall(VecResetArray(ts->vec_drdu_col)); in TSAdjointStep_Theta()
540 PetscCall(VecScale(ts->vecs_sensi[nadj], 1. / th->shift)); in TSAdjointStep_Theta()
544 if (ts->vecs_sensi2) { /* U_n */ in TSAdjointStep_Theta()
547 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); in TSAdjointStep_Theta()
549 …PetscCall(TSComputeIHessianProductFunctionUU(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_sens… in TSAdjointStep_Theta()
550 PetscCall(VecResetArray(ts->vec_sensip_col)); in TSAdjointStep_Theta()
553 …PetscCall(TSComputeIHessianProductFunctionUP(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_dir,… in TSAdjointStep_Theta()
554 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
556 PetscCall(MatMultTranspose(J, VecsDeltaLam2[nadj], ts->vecs_sensi2[nadj])); in TSAdjointStep_Theta()
557 PetscCall(VecAXPY(ts->vecs_sensi2[nadj], 1., ts->vecs_fuu[nadj])); in TSAdjointStep_Theta()
558 if (ts->vecs_fup) PetscCall(VecAXPY(ts->vecs_sensi2[nadj], 1., ts->vecs_fup[nadj])); in TSAdjointStep_Theta()
559 PetscCall(VecScale(ts->vecs_sensi2[nadj], 1. / th->shift)); in TSAdjointStep_Theta()
563 th->stage_time = ts->ptime; /* recover the old value */ in TSAdjointStep_Theta()
565 if (ts->vecs_sensip) { /* sensitivities wrt parameters */ in TSAdjointStep_Theta()
568 PetscCall(VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, ts->vec_sol)); in TSAdjointStep_Theta()
569 …PetscCall(TSComputeIJacobianP(ts, th->stage_time, ts->vec_sol, th->Xdot, -1. / (th->Theta * adjoin… in TSAdjointStep_Theta()
570 if (quadts) PetscCall(TSComputeRHSJacobianP(quadts, th->stage_time, ts->vec_sol, quadJp)); in TSAdjointStep_Theta()
571 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
572 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj])); in TSAdjointStep_Theta()
573 … PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step * th->Theta, VecsDeltaMu[nadj])); in TSAdjointStep_Theta()
576 PetscCall(VecPlaceArray(ts->vec_drdp_col, xarr)); in TSAdjointStep_Theta()
577 … PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step * th->Theta, ts->vec_drdp_col)); in TSAdjointStep_Theta()
578 PetscCall(VecResetArray(ts->vec_drdp_col)); in TSAdjointStep_Theta()
582 if (ts->vecs_sensi2p) { /* second-order */ in TSAdjointStep_Theta()
584 PetscCall(MatDenseGetColumn(ts->mat_sensip, 0, &xarr)); in TSAdjointStep_Theta()
585 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); in TSAdjointStep_Theta()
587 …scCall(TSComputeIHessianProductFunctionPU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_s… in TSAdjointStep_Theta()
588 PetscCall(VecResetArray(ts->vec_sensip_col)); in TSAdjointStep_Theta()
589 PetscCall(MatDenseRestoreColumn(ts->mat_sensip, &xarr)); in TSAdjointStep_Theta()
592 …PetscCall(TSComputeIHessianProductFunctionPP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->ve… in TSAdjointStep_Theta()
593 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
595 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam2[nadj], VecsDeltaMu2[nadj])); in TSAdjointStep_Theta()
596 … PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * th->Theta, VecsDeltaMu2[nadj])); in TSAdjointStep_Theta()
597 …if (ts->vecs_fpu) PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * th->Theta, ts->ve… in TSAdjointStep_Theta()
598 …if (ts->vecs_fpp) PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * th->Theta, ts->ve… in TSAdjointStep_Theta()
604 …PetscCall(TSComputeIJacobianP(ts, adjoint_ptime, th->X0, th->Xdot, 1. / ((th->Theta - 1.0) * adjoi… in TSAdjointStep_Theta()
606 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
607 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj])); in TSAdjointStep_Theta()
608 …PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step * (1.0 - th->Theta), VecsDeltaMu[nadj]… in TSAdjointStep_Theta()
611 PetscCall(VecPlaceArray(ts->vec_drdp_col, xarr)); in TSAdjointStep_Theta()
612 …PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step * (1.0 - th->Theta), ts->vec_drdp_col)); in TSAdjointStep_Theta()
613 PetscCall(VecResetArray(ts->vec_drdp_col)); in TSAdjointStep_Theta()
616 if (ts->vecs_sensi2p) { /* second-order */ in TSAdjointStep_Theta()
619 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); in TSAdjointStep_Theta()
621 …PetscCall(TSComputeIHessianProductFunctionPU(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_sens… in TSAdjointStep_Theta()
622 PetscCall(VecResetArray(ts->vec_sensip_col)); in TSAdjointStep_Theta()
625 …PetscCall(TSComputeIHessianProductFunctionPP(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_dir,… in TSAdjointStep_Theta()
626 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
628 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam2[nadj], VecsDeltaMu2[nadj])); in TSAdjointStep_Theta()
629 …PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * (1.0 - th->Theta), VecsDeltaMu2[nad… in TSAdjointStep_Theta()
630 …if (ts->vecs_fpu) PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * (1.0 - th->Theta)… in TSAdjointStep_Theta()
631 …if (ts->vecs_fpp) PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * (1.0 - th->Theta)… in TSAdjointStep_Theta()
638 PetscCall(TSComputeSNESJacobian(ts, th->X, J, Jpre)); /* get -f_y */ in TSAdjointStep_Theta()
641 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
643 PetscCall(VecAXPY(ts->vecs_sensi[nadj], -adjoint_time_step, VecsSensiTemp[nadj])); in TSAdjointStep_Theta()
646 PetscCall(VecPlaceArray(ts->vec_drdu_col, xarr)); in TSAdjointStep_Theta()
647 PetscCall(VecAXPY(ts->vecs_sensi[nadj], adjoint_time_step, ts->vec_drdu_col)); in TSAdjointStep_Theta()
648 PetscCall(VecResetArray(ts->vec_drdu_col)); in TSAdjointStep_Theta()
652 if (ts->vecs_sensip) { in TSAdjointStep_Theta()
655 …PetscCall(TSComputeIJacobianP(ts, th->stage_time, th->X, th->Xdot, th->shift, ts->Jacp, PETSC_FALS… in TSAdjointStep_Theta()
657 for (nadj = 0; nadj < ts->numcost; nadj++) { in TSAdjointStep_Theta()
658 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj])); in TSAdjointStep_Theta()
659 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsDeltaMu[nadj])); in TSAdjointStep_Theta()
662 PetscCall(VecPlaceArray(ts->vec_drdp_col, xarr)); in TSAdjointStep_Theta()
663 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, ts->vec_drdp_col)); in TSAdjointStep_Theta()
664 PetscCall(VecResetArray(ts->vec_drdp_col)); in TSAdjointStep_Theta()
675 static PetscErrorCode TSInterpolate_Theta(TS ts, PetscReal t, Vec X) in TSInterpolate_Theta() argument
677 TS_Theta *th = (TS_Theta *)ts->data; in TSInterpolate_Theta()
678 PetscReal dt = t - ts->ptime; in TSInterpolate_Theta()
681 PetscCall(VecCopy(ts->vec_sol, th->X)); in TSInterpolate_Theta()
687 static PetscErrorCode TSEvaluateWLTE_Theta(TS ts, NormType wnormtype, PetscInt *order, PetscReal *w… in TSEvaluateWLTE_Theta() argument
689 TS_Theta *th = (TS_Theta *)ts->data; in TSEvaluateWLTE_Theta()
690 Vec X = ts->vec_sol; /* X = solution */ in TSEvaluateWLTE_Theta()
700 if (ts->steprestart) { in TSEvaluateWLTE_Theta()
706 PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev; in TSEvaluateWLTE_Theta()
719 PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter)); in TSEvaluateWLTE_Theta()
725 static PetscErrorCode TSRollBack_Theta(TS ts) in TSRollBack_Theta() argument
727 TS_Theta *th = (TS_Theta *)ts->data; in TSRollBack_Theta()
728 TS quadts = ts->quadraturets; in TSRollBack_Theta()
731 if (quadts && ts->costintegralfwd) PetscCall(VecCopy(th->VecCostIntegral0, quadts->vec_sol)); in TSRollBack_Theta()
733 if (ts->mat_sensip) PetscCall(MatCopy(th->MatFwdSensip0, ts->mat_sensip, SAME_NONZERO_PATTERN)); in TSRollBack_Theta()
738 static PetscErrorCode TSForwardStep_Theta(TS ts) in TSForwardStep_Theta() argument
740 TS_Theta *th = (TS_Theta *)ts->data; in TSForwardStep_Theta()
741 TS quadts = ts->quadraturets; in TSForwardStep_Theta()
752 PetscCall(MatCopy(ts->mat_sensip, th->MatFwdSensip0, SAME_NONZERO_PATTERN)); in TSForwardStep_Theta()
755 PetscCall(SNESGetKSP(ts->snes, &ksp)); in TSForwardStep_Theta()
756 PetscCall(TSGetIJacobian(ts, &J, &Jpre, NULL, NULL)); in TSForwardStep_Theta()
765 … PetscCall(TSComputeIJacobian(ts, th->ptime0, th->X0, th->Xdot, th->shift, J, Jpre, PETSC_FALSE)); in TSForwardStep_Theta()
766 PetscCall(MatMatMult(J, ts->mat_sensip, MAT_REUSE_MATRIX, PETSC_DETERMINE, &MatDeltaFwdSensip)); in TSForwardStep_Theta()
770 if (ts->Jacp) { in TSForwardStep_Theta()
772 …PetscCall(TSComputeIJacobianP(ts, th->ptime0, th->X0, th->Xdot, th->shift, ts->Jacp, PETSC_FALSE)); in TSForwardStep_Theta()
773 …PetscCall(MatAXPY(MatDeltaFwdSensip, (th->Theta - 1.) / th->Theta, ts->Jacp, SUBSET_NONZERO_PATTER… in TSForwardStep_Theta()
775 PetscCall(VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, ts->vec_sol)); in TSForwardStep_Theta()
776 …PetscCall(TSComputeIJacobianP(ts, th->stage_time, ts->vec_sol, th->Xdot, th->shift, ts->Jacp, PETS… in TSForwardStep_Theta()
777 PetscCall(MatAXPY(MatDeltaFwdSensip, -1., ts->Jacp, SUBSET_NONZERO_PATTERN)); in TSForwardStep_Theta()
781 …PetscCall(TSComputeIJacobian(ts, th->stage_time, th->X, th->Xdot, th->shift, J, Jpre, PETSC_FALSE)… in TSForwardStep_Theta()
782 PetscCall(MatMatMult(J, ts->mat_sensip, MAT_REUSE_MATRIX, PETSC_DETERMINE, &MatDeltaFwdSensip)); in TSForwardStep_Theta()
786 if (ts->Jacp) { in TSForwardStep_Theta()
789 …PetscCall(TSComputeIJacobianP(ts, th->stage_time, th->X, th->Xdot, th->shift, ts->Jacp, PETSC_FALS… in TSForwardStep_Theta()
790 PetscCall(MatAXPY(MatDeltaFwdSensip, -1., ts->Jacp, SUBSET_NONZERO_PATTERN)); in TSForwardStep_Theta()
797 …PetscCall(TSComputeIJacobian(ts, th->stage_time, ts->vec_sol, th->Xdot, th->shift, J, Jpre, PETSC_… in TSForwardStep_Theta()
799 …PetscCall(TSComputeIJacobian(ts, th->stage_time, th->X, th->Xdot, th->shift, J, Jpre, PETSC_FALSE)… in TSForwardStep_Theta()
812 …PetscCall(MatTransposeMatMult(ts->mat_sensip, quadJ, MAT_REUSE_MATRIX, PETSC_DETERMINE, &th->MatIn… in TSForwardStep_Theta()
824 PetscCall(MatDenseGetColumn(ts->mat_sensip, ntlm, &xarr)); in TSForwardStep_Theta()
825 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); in TSForwardStep_Theta()
826 PetscCall(KSPSolve(ksp, VecDeltaFwdSensipCol, ts->vec_sensip_col)); in TSForwardStep_Theta()
827 PetscCall(VecResetArray(ts->vec_sensip_col)); in TSForwardStep_Theta()
828 PetscCall(MatDenseRestoreColumn(ts->mat_sensip, &xarr)); in TSForwardStep_Theta()
834 ts->reason = TSFORWARD_DIVERGED_LINEAR_SOLVE; in TSForwardStep_Theta()
835 …ll(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th tangent linear solve, linear solve f… in TSForwardStep_Theta()
847 …PetscCall(MatAXPY(ts->mat_sensip, 1, MatDeltaFwdSensip, SAME_NONZERO_PATTERN)); /* stage sensitivi… in TSForwardStep_Theta()
850 …PetscCall(MatTransposeMatMult(ts->mat_sensip, quadJ, MAT_REUSE_MATRIX, PETSC_DETERMINE, &th->MatIn… in TSForwardStep_Theta()
853 …PetscCall(MatAXPY(ts->mat_sensip, (1. - th->Theta) / th->Theta, MatDeltaFwdSensip, SAME_NONZERO_PA… in TSForwardStep_Theta()
855 PetscCall(TSComputeRHSJacobian(quadts, th->stage_time, ts->vec_sol, quadJ, NULL)); in TSForwardStep_Theta()
856 PetscCall(TSComputeRHSJacobianP(quadts, th->stage_time, ts->vec_sol, quadJp)); in TSForwardStep_Theta()
857 …PetscCall(MatTransposeMatMult(ts->mat_sensip, quadJ, MAT_REUSE_MATRIX, PETSC_DETERMINE, &th->MatIn… in TSForwardStep_Theta()
862 …if (!th->endpoint) PetscCall(MatAXPY(ts->mat_sensip, 1. / th->Theta, MatDeltaFwdSensip, SAME_NONZE… in TSForwardStep_Theta()
867 static PetscErrorCode TSForwardGetStages_Theta(TS ts, PetscInt *ns, Mat *stagesensip[]) in TSForwardGetStages_Theta() argument
869 TS_Theta *th = (TS_Theta *)ts->data; in TSForwardGetStages_Theta()
881 th->MatFwdStages[1] = ts->mat_sensip; /* stiffly accurate */ in TSForwardGetStages_Theta()
888 static PetscErrorCode TSReset_Theta(TS ts) in TSReset_Theta() argument
890 TS_Theta *th = (TS_Theta *)ts->data; in TSReset_Theta()
905 static PetscErrorCode TSAdjointReset_Theta(TS ts) in TSAdjointReset_Theta() argument
907 TS_Theta *th = (TS_Theta *)ts->data; in TSAdjointReset_Theta()
910 PetscCall(VecDestroyVecs(ts->numcost, &th->VecsDeltaLam)); in TSAdjointReset_Theta()
911 PetscCall(VecDestroyVecs(ts->numcost, &th->VecsDeltaMu)); in TSAdjointReset_Theta()
912 PetscCall(VecDestroyVecs(ts->numcost, &th->VecsDeltaLam2)); in TSAdjointReset_Theta()
913 PetscCall(VecDestroyVecs(ts->numcost, &th->VecsDeltaMu2)); in TSAdjointReset_Theta()
914 PetscCall(VecDestroyVecs(ts->numcost, &th->VecsSensiTemp)); in TSAdjointReset_Theta()
915 PetscCall(VecDestroyVecs(ts->numcost, &th->VecsSensi2Temp)); in TSAdjointReset_Theta()
919 static PetscErrorCode TSDestroy_Theta(TS ts) in TSDestroy_Theta() argument
922 PetscCall(TSReset_Theta(ts)); in TSDestroy_Theta()
923 if (ts->dm) { in TSDestroy_Theta()
924 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSTheta, DMRestrictHook_TSTheta, ts)); in TSDestroy_Theta()
925 …PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSTheta, DMSubDomainRestrictHook_TSTheta, … in TSDestroy_Theta()
927 PetscCall(PetscFree(ts->data)); in TSDestroy_Theta()
928 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetTheta_C", NULL)); in TSDestroy_Theta()
929 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetTheta_C", NULL)); in TSDestroy_Theta()
930 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetEndpoint_C", NULL)); in TSDestroy_Theta()
931 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetEndpoint_C", NULL)); in TSDestroy_Theta()
943 static PetscErrorCode SNESTSFormFunction_Theta(SNES snes, Vec x, Vec y, TS ts) in SNESTSFormFunction_Theta() argument
945 TS_Theta *th = (TS_Theta *)ts->data; in SNESTSFormFunction_Theta()
953 PetscCall(TSThetaGetX0AndXdot(ts, dm, &X0, &Xdot)); in SNESTSFormFunction_Theta()
960 dmsave = ts->dm; in SNESTSFormFunction_Theta()
961 ts->dm = dm; in SNESTSFormFunction_Theta()
962 PetscCall(TSComputeIFunction(ts, th->stage_time, x, Xdot, y, PETSC_FALSE)); in SNESTSFormFunction_Theta()
963 ts->dm = dmsave; in SNESTSFormFunction_Theta()
964 PetscCall(TSThetaRestoreX0AndXdot(ts, dm, &X0, &Xdot)); in SNESTSFormFunction_Theta()
968 static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes, Vec x, Mat A, Mat B, TS ts) in SNESTSFormJacobian_Theta() argument
970 TS_Theta *th = (TS_Theta *)ts->data; in SNESTSFormJacobian_Theta()
978 PetscCall(TSThetaGetX0AndXdot(ts, dm, NULL, &Xdot)); in SNESTSFormJacobian_Theta()
980 dmsave = ts->dm; in SNESTSFormJacobian_Theta()
981 ts->dm = dm; in SNESTSFormJacobian_Theta()
982 PetscCall(TSComputeIJacobian(ts, th->stage_time, x, Xdot, shift, A, B, PETSC_FALSE)); in SNESTSFormJacobian_Theta()
983 ts->dm = dmsave; in SNESTSFormJacobian_Theta()
984 PetscCall(TSThetaRestoreX0AndXdot(ts, dm, NULL, &Xdot)); in SNESTSFormJacobian_Theta()
988 static PetscErrorCode TSForwardSetUp_Theta(TS ts) in TSForwardSetUp_Theta() argument
990 TS_Theta *th = (TS_Theta *)ts->data; in TSForwardSetUp_Theta()
991 TS quadts = ts->quadraturets; in TSForwardSetUp_Theta()
995 th->num_tlm = ts->num_parameters; in TSForwardSetUp_Theta()
996 PetscCall(MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &th->MatDeltaFwdSensip)); in TSForwardSetUp_Theta()
1002 PetscCall(MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &th->MatFwdSensip0)); in TSForwardSetUp_Theta()
1004 PetscCall(VecDuplicate(ts->vec_sol, &th->VecDeltaFwdSensipCol)); in TSForwardSetUp_Theta()
1008 static PetscErrorCode TSForwardReset_Theta(TS ts) in TSForwardReset_Theta() argument
1010 TS_Theta *th = (TS_Theta *)ts->data; in TSForwardReset_Theta()
1011 TS quadts = ts->quadraturets; in TSForwardReset_Theta()
1024 static PetscErrorCode TSSetUp_Theta(TS ts) in TSSetUp_Theta() argument
1026 TS_Theta *th = (TS_Theta *)ts->data; in TSSetUp_Theta()
1027 TS quadts = ts->quadraturets; in TSSetUp_Theta()
1031 if (!th->VecCostIntegral0 && quadts && ts->costintegralfwd) { /* back up cost integral */ in TSSetUp_Theta()
1034 if (!th->X) PetscCall(VecDuplicate(ts->vec_sol, &th->X)); in TSSetUp_Theta()
1035 if (!th->Xdot) PetscCall(VecDuplicate(ts->vec_sol, &th->Xdot)); in TSSetUp_Theta()
1036 if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0)); in TSSetUp_Theta()
1037 if (th->endpoint) PetscCall(VecDuplicate(ts->vec_sol, &th->affine)); in TSSetUp_Theta()
1040 th->shift = 1 / (th->Theta * ts->time_step); in TSSetUp_Theta()
1042 PetscCall(TSGetDM(ts, &ts->dm)); in TSSetUp_Theta()
1043 PetscCall(DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSTheta, DMRestrictHook_TSTheta, ts)); in TSSetUp_Theta()
1044 …PetscCall(DMSubDomainHookAdd(ts->dm, DMSubDomainHook_TSTheta, DMSubDomainRestrictHook_TSTheta, ts)… in TSSetUp_Theta()
1046 PetscCall(TSGetAdapt(ts, &ts->adapt)); in TSSetUp_Theta()
1047 PetscCall(TSAdaptCandidatesClear(ts->adapt)); in TSSetUp_Theta()
1048 PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match)); in TSSetUp_Theta()
1050 if (!th->vec_sol_prev) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev)); in TSSetUp_Theta()
1051 if (!th->vec_lte_work) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work)); in TSSetUp_Theta()
1053 PetscCall(TSGetSNES(ts, &ts->snes)); in TSSetUp_Theta()
1055 ts->stifflyaccurate = (!th->endpoint && th->Theta != 1.0) ? PETSC_FALSE : PETSC_TRUE; in TSSetUp_Theta()
1059 static PetscErrorCode TSAdjointSetUp_Theta(TS ts) in TSAdjointSetUp_Theta() argument
1061 TS_Theta *th = (TS_Theta *)ts->data; in TSAdjointSetUp_Theta()
1064 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &th->VecsDeltaLam)); in TSAdjointSetUp_Theta()
1065 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &th->VecsSensiTemp)); in TSAdjointSetUp_Theta()
1066 …if (ts->vecs_sensip) PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &th->VecsDeltaMu)… in TSAdjointSetUp_Theta()
1067 if (ts->vecs_sensi2) { in TSAdjointSetUp_Theta()
1068 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &th->VecsDeltaLam2)); in TSAdjointSetUp_Theta()
1069 PetscCall(VecDuplicateVecs(ts->vecs_sensi2[0], ts->numcost, &th->VecsSensi2Temp)); in TSAdjointSetUp_Theta()
1071 if (!ts->ihessianproduct_fuu) ts->vecs_fuu = ts->vecs_guu; in TSAdjointSetUp_Theta()
1072 if (!ts->ihessianproduct_fup) ts->vecs_fup = ts->vecs_gup; in TSAdjointSetUp_Theta()
1074 if (ts->vecs_sensi2p) { in TSAdjointSetUp_Theta()
1075 PetscCall(VecDuplicateVecs(ts->vecs_sensi2p[0], ts->numcost, &th->VecsDeltaMu2)); in TSAdjointSetUp_Theta()
1077 if (!ts->ihessianproduct_fpu) ts->vecs_fpu = ts->vecs_gpu; in TSAdjointSetUp_Theta()
1078 if (!ts->ihessianproduct_fpp) ts->vecs_fpp = ts->vecs_gpp; in TSAdjointSetUp_Theta()
1083 static PetscErrorCode TSSetFromOptions_Theta(TS ts, PetscOptionItems PetscOptionsObject) in TSSetFromOptions_Theta() argument
1085 TS_Theta *th = (TS_Theta *)ts->data; in TSSetFromOptions_Theta()
1098 static PetscErrorCode TSView_Theta(TS ts, PetscViewer viewer) in TSView_Theta() argument
1100 TS_Theta *th = (TS_Theta *)ts->data; in TSView_Theta()
1112 static PetscErrorCode TSThetaGetTheta_Theta(TS ts, PetscReal *theta) in TSThetaGetTheta_Theta() argument
1114 TS_Theta *th = (TS_Theta *)ts->data; in TSThetaGetTheta_Theta()
1121 static PetscErrorCode TSThetaSetTheta_Theta(TS ts, PetscReal theta) in TSThetaSetTheta_Theta() argument
1123 TS_Theta *th = (TS_Theta *)ts->data; in TSThetaSetTheta_Theta()
1126 …PetscCheck(theta > 0 && theta <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "T… in TSThetaSetTheta_Theta()
1132 static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts, PetscBool *endpoint) in TSThetaGetEndpoint_Theta() argument
1134 TS_Theta *th = (TS_Theta *)ts->data; in TSThetaGetEndpoint_Theta()
1141 static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts, PetscBool flg) in TSThetaSetEndpoint_Theta() argument
1143 TS_Theta *th = (TS_Theta *)ts->data; in TSThetaSetEndpoint_Theta()
1151 static PetscErrorCode TSComputeLinearStability_Theta(TS ts, PetscReal xr, PetscReal xi, PetscReal *… in TSComputeLinearStability_Theta() argument
1154 TS_Theta *th = (TS_Theta *)ts->data; in TSComputeLinearStability_Theta()
1164 static PetscErrorCode TSGetStages_Theta(TS ts, PetscInt *ns, Vec *Y[]) in TSGetStages_Theta() argument
1166 TS_Theta *th = (TS_Theta *)ts->data; in TSGetStages_Theta()
1178 th->Stages[1] = ts->vec_sol; /* stiffly accurate */ in TSGetStages_Theta()
1233 PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts) in TSCreate_Theta() argument
1238 ts->ops->reset = TSReset_Theta; in TSCreate_Theta()
1239 ts->ops->adjointreset = TSAdjointReset_Theta; in TSCreate_Theta()
1240 ts->ops->destroy = TSDestroy_Theta; in TSCreate_Theta()
1241 ts->ops->view = TSView_Theta; in TSCreate_Theta()
1242 ts->ops->setup = TSSetUp_Theta; in TSCreate_Theta()
1243 ts->ops->adjointsetup = TSAdjointSetUp_Theta; in TSCreate_Theta()
1244 ts->ops->adjointreset = TSAdjointReset_Theta; in TSCreate_Theta()
1245 ts->ops->step = TSStep_Theta; in TSCreate_Theta()
1246 ts->ops->interpolate = TSInterpolate_Theta; in TSCreate_Theta()
1247 ts->ops->evaluatewlte = TSEvaluateWLTE_Theta; in TSCreate_Theta()
1248 ts->ops->rollback = TSRollBack_Theta; in TSCreate_Theta()
1249 ts->ops->resizeregister = TSResizeRegister_Theta; in TSCreate_Theta()
1250 ts->ops->setfromoptions = TSSetFromOptions_Theta; in TSCreate_Theta()
1251 ts->ops->snesfunction = SNESTSFormFunction_Theta; in TSCreate_Theta()
1252 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; in TSCreate_Theta()
1254 ts->ops->linearstability = TSComputeLinearStability_Theta; in TSCreate_Theta()
1256 ts->ops->getstages = TSGetStages_Theta; in TSCreate_Theta()
1257 ts->ops->adjointstep = TSAdjointStep_Theta; in TSCreate_Theta()
1258 ts->ops->adjointintegral = TSAdjointCostIntegral_Theta; in TSCreate_Theta()
1259 ts->ops->forwardintegral = TSForwardCostIntegral_Theta; in TSCreate_Theta()
1260 ts->default_adapt_type = TSADAPTNONE; in TSCreate_Theta()
1262 ts->ops->forwardsetup = TSForwardSetUp_Theta; in TSCreate_Theta()
1263 ts->ops->forwardreset = TSForwardReset_Theta; in TSCreate_Theta()
1264 ts->ops->forwardstep = TSForwardStep_Theta; in TSCreate_Theta()
1265 ts->ops->forwardgetstages = TSForwardGetStages_Theta; in TSCreate_Theta()
1267 ts->usessnes = PETSC_TRUE; in TSCreate_Theta()
1270 ts->data = (void *)th; in TSCreate_Theta()
1280 …PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetTheta_C", TSThetaGetTheta_Theta)); in TSCreate_Theta()
1281 …PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetTheta_C", TSThetaSetTheta_Theta)); in TSCreate_Theta()
1282 …PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetEndpoint_C", TSThetaGetEndpoint_T… in TSCreate_Theta()
1283 …PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetEndpoint_C", TSThetaSetEndpoint_T… in TSCreate_Theta()
1305 PetscErrorCode TSThetaGetTheta(TS ts, PetscReal *theta) in TSThetaGetTheta() argument
1308 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); in TSThetaGetTheta()
1310 PetscUseMethod(ts, "TSThetaGetTheta_C", (TS, PetscReal *), (ts, theta)); in TSThetaGetTheta()
1330 PetscErrorCode TSThetaSetTheta(TS ts, PetscReal theta) in TSThetaSetTheta() argument
1333 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); in TSThetaSetTheta()
1334 PetscTryMethod(ts, "TSThetaSetTheta_C", (TS, PetscReal), (ts, theta)); in TSThetaSetTheta()
1353 PetscErrorCode TSThetaGetEndpoint(TS ts, PetscBool *endpoint) in TSThetaGetEndpoint() argument
1356 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); in TSThetaGetEndpoint()
1358 PetscUseMethod(ts, "TSThetaGetEndpoint_C", (TS, PetscBool *), (ts, endpoint)); in TSThetaGetEndpoint()
1378 PetscErrorCode TSThetaSetEndpoint(TS ts, PetscBool flg) in TSThetaSetEndpoint() argument
1381 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); in TSThetaSetEndpoint()
1382 PetscTryMethod(ts, "TSThetaSetEndpoint_C", (TS, PetscBool), (ts, flg)); in TSThetaSetEndpoint()
1391 static PetscErrorCode TSSetUp_BEuler(TS ts) in TSSetUp_BEuler() argument
1393 TS_Theta *th = (TS_Theta *)ts->data; in TSSetUp_BEuler()
1396 …PetscCheck(th->Theta == 1.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_OPT_OVERWRITE, "Can not c… in TSSetUp_BEuler()
1397 …PetscCheck(!th->endpoint, PetscObjectComm((PetscObject)ts), PETSC_ERR_OPT_OVERWRITE, "Can not chan… in TSSetUp_BEuler()
1398 PetscCall(TSSetUp_Theta(ts)); in TSSetUp_BEuler()
1402 static PetscErrorCode TSView_BEuler(TS ts, PetscViewer viewer) in TSView_BEuler() argument
1418 PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts) in TSCreate_BEuler() argument
1421 PetscCall(TSCreate_Theta(ts)); in TSCreate_BEuler()
1422 PetscCall(TSThetaSetTheta(ts, 1.0)); in TSCreate_BEuler()
1423 PetscCall(TSThetaSetEndpoint(ts, PETSC_FALSE)); in TSCreate_BEuler()
1424 ts->ops->setup = TSSetUp_BEuler; in TSCreate_BEuler()
1425 ts->ops->view = TSView_BEuler; in TSCreate_BEuler()
1429 static PetscErrorCode TSSetUp_CN(TS ts) in TSSetUp_CN() argument
1431 TS_Theta *th = (TS_Theta *)ts->data; in TSSetUp_CN()
1434 …PetscCheck(th->Theta == 0.5, PetscObjectComm((PetscObject)ts), PETSC_ERR_OPT_OVERWRITE, "Can not c… in TSSetUp_CN()
1435 …PetscCheck(th->endpoint, PetscObjectComm((PetscObject)ts), PETSC_ERR_OPT_OVERWRITE, "Can not chang… in TSSetUp_CN()
1436 PetscCall(TSSetUp_Theta(ts)); in TSSetUp_CN()
1440 static PetscErrorCode TSView_CN(TS ts, PetscViewer viewer) in TSView_CN() argument
1461 PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts) in TSCreate_CN() argument
1464 PetscCall(TSCreate_Theta(ts)); in TSCreate_CN()
1465 PetscCall(TSThetaSetTheta(ts, 0.5)); in TSCreate_CN()
1466 PetscCall(TSThetaSetEndpoint(ts, PETSC_TRUE)); in TSCreate_CN()
1467 ts->ops->setup = TSSetUp_CN; in TSCreate_CN()
1468 ts->ops->view = TSView_CN; in TSCreate_CN()