On Jul 26, 2017, at 8:34 AM, Stefano Zampini <[email protected]<mailto:[email protected]>> wrote:
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 This is the reason why we use ts->total_steps (global index) instead of ts->steps (local index). There are use cases that a forward simulation consists of multiple consecutive TSSolve(), e.g. power_grid/stability_9bus/ex9busadj.c And ts->total_steps was designed to support these cases. It allows us to compute the adjoint by doing TSSolve() // forward propagation for ts->total_steps=[0, n_steps_1] TSSolve() // forward propagation for ts->total_steps=[n_steps_1, n_steps_1+n_steps_2] TSAdjointSolve() // backward propagation for ts->total_steps=[n_steps_1 + n_steps_2, 0] The output of TSMonitor/TSAdjointMonitor can reflect this workflow nicely. For each TSSolve, ts->steps always starts from 0, if we use it for trajectory, the two TSSolve() would cause conflicting step numbers. Hong 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 <[email protected]<mailto:[email protected]>>: 2017-07-25 22:45 GMT+03:00 Zhang, Hong <[email protected]<mailto:[email protected]>>: > > 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
