this is regarding ts->total_steps. What is the actual purpose of it? You claim it is a unique identifier, but actually TSAdjointStep decrements it. We should come to an agreement on it, as it is used to monitoring a TSSolve call, and to save the trajectory, as per commit
https://bitbucket.org/petsc/petsc/commits/938c94ca59bdacf9e5030aef0fd03ca8c4ad6971 Right now, if we do TSSolve() // with number of steps n_steps_1 TSSolve() // with number of steps n_steps_2 the output of a monitor will be 0 ... .. n_steps_1 // here starts the monitoring of the second solve n_steps_1 .... n_steps_1 + n_steps_2 Instead, if we do TSSolve() TSAdjointSolve() TSSolve() the output will be 0 ... .. n_steps_1 // here starts the monitoring of the adjoint solve n_steps_1 ... 0 // here starts the monitoring of the second solve 0 .... n_steps_1 Why not using ts->steps instead of ts->total_steps in TSMonitor and TSSetTrajectory? We should be really careful when modifying interface code. 2017-07-26 13:11 GMT+03:00 Stefano Zampini <stefano.zamp...@gmail.com>: > > > 2017-07-25 22:45 GMT+03:00 Zhang, Hong <hongzh...@anl.gov>: > >> >> >> > >> > Actually, it is not usable: take for example TSTrajectoryGet() with a >> given step number: then, if I want to use BASIC it will always return the >> last step recorded from the previous TSSolve (this behavior is common to >> all the TSTrajectory implementantions) >> >> This is not the expected behavior. Can you provide more details on how >> you used TSTrajectoryGet()? How did you know it is the last step recorded >> from the previous TSSolve? >> >> > Hong, > > it's clear from the code, since in every implementation of > TSTrajectoryGet_XXX TSGetTotalSteps is overwriting the stepnumber passed in > by the user > A patch is attached, that removes this extra call. Could you please check > that the patch does not break the adjoint examples? > The patch should be fine, since TSTrajectoryGet is just used in > TSAdjointSolve and it is always called using ts->total_steps as stepnum. > > [szampini@localhost petsc]$ git grep "TSTrajectoryGet(" > include/petscts.h:PETSC_EXTERN PetscErrorCode TSTrajectoryGet(TSTrajectory, > TS,PetscInt,PetscReal*); > src/ts/interface/ts.c: ierr = TSTrajectoryGet(ts-> > trajectory,ts,ts->total_steps,&ts->ptime);CHKERRQ(ierr); > src/ts/interface/ts.c: ierr = TSTrajectoryGet(ts-> > trajectory,ts,ts->total_steps,&ts->ptime);CHKERRQ(ierr); > src/ts/trajectory/interface/traj.c:PetscErrorCode > TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time) > > >> >> total_steps is automatically incremented in TSSolve(), which serves as a >> unique identifier for a time step. TSTrajectory is marked as 'developer >> level'. The implementation is quite tricky. I am willing to help if you >> have some specific needs. >> >> > The continuous adjoint code that I have uses two different TS solvers, one > for forward, one for backward. To compute the gradient, I need to inquire > the trajectory of the forward solver in a post step routine and compute the > gradient. Inquiring the trajectory by the time argument and not only by the > stepnum would be the best, so that one can (in principle) use adaptive time > stepping in both solvers, without worrying about the matching time steps. > For this, it would be great if we can have a TSTrajectory that computes on > the fly if the requested time is not present in the recorded history. > For now, I can write my own map from forward time to forward step number > and inquire the trajectory just with the step number. > > > >> Thanks, >> Hong (Mr.) >> >> > >> > I would really like to reuse TSTrajectory for a continuous adjoint >> problem I'm currently working on, and I believe this class should be >> somewhat refactored. Here are a couple of examples from the public API >> > >> > PetscErrorCode TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt >> stepnum,PetscReal *time) >> > { >> > PetscErrorCode ierr; >> > >> > PetscFunctionBegin; >> > if (!tj) >> > SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TS >> solver did not save trajectory"); >> > >> > .... >> > >> > PetscErrorCode TSTrajectorySet(TSTrajectory tj,TS ts,PetscInt >> stepnum,PetscReal time,Vec X) >> > { >> > PetscErrorCode ierr; >> > >> > PetscFunctionBegin; >> > if (!tj) PetscFunctionReturn(0); >> > >> > We can use this thread to start a discussion on how it can be >> refactored, and I'm willing to contribute. >> > >> > I personally believe TSTrajectory should not be exposed in the public >> API (at least put in the private headers) in the upcoming release. >> > >> > Thanks, >> > -- >> > Stefano >> >> > > > -- > Stefano > -- Stefano
