<div dir="ltr">I don't have a real API, just a sketch of a possible general framework (with holes....)<div><br></div><div>First, I prefer having TSCreateAdjointTS(ts,&ats) (and possibly SNESGetAdjointSNES for optimization with nonlinear time-independent problems) that allows users to replace their Jacobian evaluation routines.<div><br></div><div>Then, there's the need of carrying over the design vector to TS (and SNES), and this can be accomplished by a DM<br><div><br></div><div>// User code</div><div>DMGetDMDesign(dm,&ddm)</div><div>DMSetApplicationContext(ddm,userctx); // attach user context</div><div>DMDesignSetSetUp(ddm,(PetscErrorCode*)(DM,Vec)) //setup user context<br></div><div>DMDesignSetGlobalToLocal //maybe other callbacks ???</div><div><div>DMDesignSetDesignVec(ddm,currentdesign) //current design vector, can be updated by higher level routines<br></div></div><div><br></div><div>Then, within PETSc, we can do</div><div><br></div><div>TSGetDM(ts,&dm)</div><div>DMGetDMTS(dm,&tdm)<br></div><div>DMGetDMDesign(dm,&ddm)</div><div>DMTSSetDMDesign(tdm,ddm) // without taking reference on ddm, just passing callbacks from DMDesign to private DMTS ops </div><div><br></div><div>and for gradient based optimization and sensitivity for DAEs, we need to implement the equivalents of TSSet{Gradient|Hessian}{DAE|IC} in DMTS. DMDesign should provide an API with callbacks for computing gradients (and possibly hessian terms) of residual evaluations wrt the parameters.</div><div><br></div><div>Similarly for a parameter dependent SNES</div><div><div><br></div><div>SNESGetDM(snes,&dm)</div><div>DMGetDMSNES(dm,&kdm)<br></div><div>DMGetDMDesign(dm,&ddm)</div><div>DMSNESSetDMDesign(kdm,ddm) // without taking reference on ddm, just passing callbacks from DMDesign to private DMSNES ops</div></div><div><br></div><div>Then, we need a consistent way of expressing objective functions (DMTAOObjective???), and embed quadratures inside TSSolve(). Thoughts?</div><div><br></div><div><br></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-10-17 22:02 GMT+03:00 Jed Brown <span dir="ltr"><<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="HOEnZb"><div class="h5">Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> writes:<br>
<br>
>> On Oct 17, 2017, at 10:41 AM, Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>> wrote:<br>
>><br>
>> Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> writes:<br>
>><br>
>>>>>> My<br>
>>>>>> preference is for TSCreateAdjointTS() which the user can customize to<br>
>>>>>> use arbitrary time integration methods (thus becoming a<br>
>>>>>> "continuous-in-time" adjoint) and to use different spatial<br>
>>>>>> discretizations (for consistency with the adjoint PDE).<br>
>>>>><br>
>>>>> I understand, I think TSCreateAdjointTS() is not the right model. I think that the correct model is, for example<br>
>>>>><br>
>>>>> PetscSensiCreate()<br>
>>>>> PetscSensiSetObjective.... or whatever etc etc<br>
>>>>> PetscSensiGetIntegrator(,&TS)<br>
>>>>> PetscSensiGetAdjointIntegrator<wbr>(,&TS).<br>
>>>><br>
>>>> How would PetscSensiGetAdjointIntegrator be implemented? Since it needs<br>
>>>> to be able to create a discrete adjoint (which depends on the forward TS<br>
>>>> implementation), the control flow *must* go through TS. What does that<br>
>>>> interface look like?<br>
>>><br>
>>> I am not sure what you mean.<br>
>>><br>
>>> PetscSensi will contain the information about objective function etc as well as a pointer to the trajectory history and (of course) a pointer to the original TS. All the auxiliary vectors etc related to sensitivities.<br>
>>><br>
>>> The adjoint TS that is returned has the calls need to get access to the forward solution that it needs etc from the PetscSensi object (which likely will call the forward TS for missing timestep values etc. But the calls go through the PetscSensi object from the adjoint TS, they don't go directly to the forward TS.)<br>
>><br>
>> To be concrete with TSRK, PetscSensi does not have access to the Butcher<br>
>> table for the forward method because that only exists in<br>
>> src/ts/impls/explicit/rk/rk.c. Since the discrete adjoint depends on<br>
>> this detail of the forward method, there MUST be some code somewhere<br>
>> that depends on the forward Butcher table. Since that data ONLY EXISTS<br>
>> in rk.c, whatever it is that PetscSensi does, there must be some code in<br>
>> rk.c to define the adjoint method.<br>
><br>
> Of course, that would be incorporated in what is an essentially a TSRKADJOINT implementation that does the adjoint integration. Note that (1) this is simpler than the regular TSRK because it only needs to handle linear with given time steps (2) it also computes along sensitivity information with calls to TSAdjointComputeDRDPFunction() etc<br>
><br>
> Here is what the code looks like now. Just needs minor refactoring to make it match my proposal.<br>
<br>
</div></div>Right, so this function needs to call the dF/dp (RHSJacobian), dr/dy,<br>
and dr/dp functions. You're proposing that those would be set on<br>
PetscSensi, but this code below will continue to live in rk.c after<br>
refactoring (rather than in PetscSensi).<br>
<br>
Since PetscSensi is not inside TS, there needs to be a public TS<br>
interface ("developer level" if you like, but still public). I want to<br>
spec out that essential bit of public TS interface.<br>
<br>
Note that I proposed putting the sensitivity update stuff below into<br>
post-stage functions. If that is possible, the bulk of the function<br>
could actually live in PetscSensi. Or it could be packaged with TS and<br>
reusable across TS implementations. I still want to start by spec'ing<br>
out what the TS public interface needs to able to support this<br>
PetscSensi object.<br>
<div><div class="h5"><br>
> static PetscErrorCode TSAdjointStep_RK(TS ts)<br>
> {<br>
> TS_RK *rk = (TS_RK*)ts->data;<br>
> RKTableau tab = rk->tableau;<br>
> const PetscInt s = tab->s;<br>
> const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;<br>
> PetscScalar *w = rk->work;<br>
> Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;<br>
> PetscInt i,j,nadj;<br>
> PetscReal t = ts->ptime;<br>
> PetscErrorCode ierr;<br>
> PetscReal h = ts->time_step;<br>
> Mat J,Jp;<br>
><br>
> PetscFunctionBegin;<br>
> rk->status = TS_STEP_INCOMPLETE;<br>
> for (i=s-1; i>=0; i--) {<br>
> rk->stage_time = t + h*(1.0-c[i]);<br>
> for (nadj=0; nadj<ts->numcost; nadj++) {<br>
> ierr = VecCopy(ts->vecs_sensi[nadj],<wbr>VecSensiTemp[nadj]);CHKERRQ(<wbr>ierr);<br>
> ierr = VecScale(VecSensiTemp[nadj],-<wbr>h*b[i]);CHKERRQ(ierr);<br>
> for (j=i+1; j<s; j++) {<br>
> ierr = VecAXPY(VecSensiTemp[nadj],-h*<wbr>A[j*s+i],VecDeltaLam[nadj*s+j]<wbr>);CHKERRQ(ierr);<br>
> }<br>
> }<br>
> /* Stage values of lambda */<br>
> ierr = TSGetRHSJacobian(ts,&J,&Jp,<wbr>NULL,NULL);CHKERRQ(ierr);<br>
> ierr = TSComputeRHSJacobian(ts,rk-><wbr>stage_time,Y[i],J,Jp);CHKERRQ(<wbr>ierr);<br>
> if (ts->vec_costintegral) {<br>
> ierr = TSAdjointComputeDRDYFunction(<wbr>ts,rk->stage_time,Y[i],ts-><wbr>vecs_drdy);CHKERRQ(ierr);<br>
> }<br>
> for (nadj=0; nadj<ts->numcost; nadj++) {<br>
> ierr = MatMultTranspose(J,<wbr>VecSensiTemp[nadj],<wbr>VecDeltaLam[nadj*s+i]);<wbr>CHKERRQ(ierr);<br>
> if (ts->vec_costintegral) {<br>
> ierr = VecAXPY(VecDeltaLam[nadj*s+i],<wbr>-h*b[i],ts->vecs_drdy[nadj]);<wbr>CHKERRQ(ierr);<br>
> }<br>
> }<br>
><br>
> /* Stage values of mu */<br>
> if(ts->vecs_sensip) {<br>
> ierr = TSAdjointComputeRHSJacobian(<wbr>ts,rk->stage_time,Y[i],ts-><wbr>Jacp);CHKERRQ(ierr);<br>
> if (ts->vec_costintegral) {<br>
> ierr = TSAdjointComputeDRDPFunction(<wbr>ts,rk->stage_time,Y[i],ts-><wbr>vecs_drdp);CHKERRQ(ierr);<br>
> }<br>
><br>
> for (nadj=0; nadj<ts->numcost; nadj++) {<br>
> ierr = MatMultTranspose(ts->Jacp,<wbr>VecSensiTemp[nadj],VecDeltaMu[<wbr>nadj*s+i]);CHKERRQ(ierr);<br>
> if (ts->vec_costintegral) {<br>
> ierr = VecAXPY(VecDeltaMu[nadj*s+i],-<wbr>h*b[i],ts->vecs_drdp[nadj]);<wbr>CHKERRQ(ierr);<br>
> }<br>
> }<br>
> }<br>
> }<br>
><br>
> for (j=0; j<s; j++) w[j] = 1.0;<br>
> for (nadj=0; nadj<ts->numcost; nadj++) {<br>
> ierr = VecMAXPY(ts->vecs_sensi[nadj],<wbr>s,w,&VecDeltaLam[nadj*s]);<wbr>CHKERRQ(ierr);<br>
> if(ts->vecs_sensip) {<br>
> ierr = VecMAXPY(ts->vecs_sensip[nadj]<wbr>,s,w,&VecDeltaMu[nadj*s]);<wbr>CHKERRQ(ierr);<br>
> }<br>
> }<br>
> rk->status = TS_STEP_COMPLETE;<br>
> PetscFunctionReturn(0);<br>
> }<br>
><br>
>><br>
>> What is that code and what is the call stack to reach it?<br>
><br>
> Now you ask, how does the PetscSensi know to use that discrete adjoint integrator and return the right TS (instead of an empty one requiring the user to set the type)? This is simple, each TS can return<br>
> its discrete adjoint integrator when requested, sure mechanically it<br>
> is like a hidden TSGetAdjointTS() but from the public API<br>
<br>
</div></div>It's public TS API: PetscSensi will include petscts.h, not<br>
petsc/private/tsimpl.h. Let's work out what that TS API needs to be and<br>
then we can discuss whether it is "clearer" to set that stuff on a<br>
helper object. If all we have are pass-through interfaces from<br>
PetscSensi to TS, then it's probably not useful. If we have a<br>
meaningfully higher level interface, then it may be useful.<br>
<div class="HOEnZb"><div class="h5"><br>
> point of view it is much clearer to come from the PetscSensi. You tell<br>
> the PetscSensi all the information you need about the derivatives, you<br>
> don't tell the TS. You are also right that not EVERYTHING related to<br>
> sensitivities is only in the PetscSensi code since the TSRKADJOINT<br>
> object knows about sensitivities that it carries around (and is in<br>
> rk.c) but from the user point of view only the PetscSensi knows about<br>
> sensitivities.<br>
><br>
> Of course I could still be missing stuff<br>
><br>
> Barry<br>
><br>
><br>
>><br>
>>> Both proposed models shove one heck of a lot of stuff into TS; I think it would be a better design and better user API if we don't do that; we reserve the TS for being able to integrate something (and that is about all). Maybe it is not possible but I don't see why it is not possible.<br>
>>><br>
>>> Since I don't know the details I asked for the people who do know the details (Hong and Stefano) to try to develop such an API (or clearly articular why it is impossible.)<br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature">Stefano</div>
</div>