xref: /petsc/src/ts/tutorials/ex20adj.c (revision 5f80ce2ab25dff0f4601e710601cbbcecf323266)
1c4762a1bSJed Brown static char help[] = "Performs adjoint sensitivity analysis for the van der Pol equation.\n";
2c4762a1bSJed Brown 
3c4762a1bSJed Brown /*
4c4762a1bSJed Brown    Concepts: TS^time-dependent nonlinear problems
5c4762a1bSJed Brown    Concepts: TS^van der Pol equation DAE equivalent
6c4762a1bSJed Brown    Concepts: TS^adjoint sensitivity analysis
7c4762a1bSJed Brown    Processors: 1
8c4762a1bSJed Brown */
9c4762a1bSJed Brown /* ------------------------------------------------------------------------
10c4762a1bSJed Brown 
11c4762a1bSJed Brown    This program solves the van der Pol DAE ODE equivalent
12c4762a1bSJed Brown       [ u_1' ] = [          u_2                ]  (2)
13c4762a1bSJed Brown       [ u_2' ]   [ \mu ((1 - u_1^2) u_2 - u_1) ]
14c4762a1bSJed Brown    on the domain 0 <= x <= 1, with the boundary conditions
15c4762a1bSJed Brown        u_1(0) = 2, u_2(0) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2),
16c4762a1bSJed Brown    and
17c4762a1bSJed Brown        \mu = 10^6 ( y'(0) ~ -0.6666665432100101).,
18c4762a1bSJed Brown    and computes the sensitivities of the final solution w.r.t. initial conditions and parameter \mu with the implicit theta method and its discrete adjoint.
19c4762a1bSJed Brown 
20c4762a1bSJed Brown    Notes:
21c4762a1bSJed Brown    This code demonstrates the TSAdjoint interface to a DAE system.
22c4762a1bSJed Brown 
23c4762a1bSJed Brown    The user provides the implicit right-hand-side function
245ab1ac2bSHong Zhang    [ F(u',u,t) ] = [u' - f(u,t)] = [ u_1'] - [        u_2             ]
25c4762a1bSJed Brown                                    [ u_2']   [ \mu ((1-u_1^2)u_2-u_1) ]
26c4762a1bSJed Brown 
275ab1ac2bSHong Zhang    and the Jacobian of F (from the PETSc user manual)
28c4762a1bSJed Brown 
295ab1ac2bSHong Zhang               dF   dF
305ab1ac2bSHong Zhang    J(F) = a * -- + --
31c4762a1bSJed Brown               du'  du
32c4762a1bSJed Brown 
335ab1ac2bSHong Zhang    and the JacobianP of the explicit right-hand side of (2) f(u,t) ( which is equivalent to -F(0,u,t)).
34c4762a1bSJed Brown    df   [       0               ]
35c4762a1bSJed Brown    -- = [                       ]
36c4762a1bSJed Brown    dp   [ (1 - u_1^2) u_2 - u_1 ].
37c4762a1bSJed Brown 
38c4762a1bSJed Brown    See ex20.c for more details on the Jacobian.
39c4762a1bSJed Brown 
40c4762a1bSJed Brown   ------------------------------------------------------------------------- */
41c4762a1bSJed Brown #include <petscts.h>
42c4762a1bSJed Brown #include <petsctao.h>
43c4762a1bSJed Brown 
44c4762a1bSJed Brown typedef struct _n_User *User;
45c4762a1bSJed Brown struct _n_User {
46c4762a1bSJed Brown   PetscReal mu;
47c4762a1bSJed Brown   PetscReal next_output;
48c4762a1bSJed Brown 
49c4762a1bSJed Brown   /* Sensitivity analysis support */
50c4762a1bSJed Brown   PetscInt  steps;
51c4762a1bSJed Brown   PetscReal ftime;
52c4762a1bSJed Brown   Mat       A;                   /* Jacobian matrix */
53c4762a1bSJed Brown   Mat       Jacp;                /* JacobianP matrix */
54c4762a1bSJed Brown   Vec       U,lambda[2],mup[2];  /* adjoint variables */
55c4762a1bSJed Brown };
56c4762a1bSJed Brown 
57c4762a1bSJed Brown /* ----------------------- Explicit form of the ODE  -------------------- */
58c4762a1bSJed Brown 
59c4762a1bSJed Brown static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
60c4762a1bSJed Brown {
61c4762a1bSJed Brown   User              user = (User)ctx;
62c4762a1bSJed Brown   PetscScalar       *f;
63c4762a1bSJed Brown   const PetscScalar *u;
64c4762a1bSJed Brown 
65c4762a1bSJed Brown   PetscFunctionBeginUser;
66*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(U,&u));
67*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(F,&f));
68c4762a1bSJed Brown   f[0] = u[1];
69c4762a1bSJed Brown   f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
70*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(U,&u));
71*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(F,&f));
72c4762a1bSJed Brown   PetscFunctionReturn(0);
73c4762a1bSJed Brown }
74c4762a1bSJed Brown 
75c4762a1bSJed Brown static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
76c4762a1bSJed Brown {
77c4762a1bSJed Brown   User              user = (User)ctx;
78c4762a1bSJed Brown   PetscReal         mu   = user->mu;
79c4762a1bSJed Brown   PetscInt          rowcol[] = {0,1};
80c4762a1bSJed Brown   PetscScalar       J[2][2];
81c4762a1bSJed Brown   const PetscScalar *u;
82c4762a1bSJed Brown 
83c4762a1bSJed Brown   PetscFunctionBeginUser;
84*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(U,&u));
85c4762a1bSJed Brown   J[0][0] = 0;
86c4762a1bSJed Brown   J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
87c4762a1bSJed Brown   J[0][1] = 1.0;
88c4762a1bSJed Brown   J[1][1] = mu*(1.0-u[0]*u[0]);
89*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES));
90*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
91*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
92c4762a1bSJed Brown   if (A != B) {
93*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
94*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
95c4762a1bSJed Brown   }
96*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(U,&u));
97c4762a1bSJed Brown   PetscFunctionReturn(0);
98c4762a1bSJed Brown }
99c4762a1bSJed Brown 
100c4762a1bSJed Brown /* ----------------------- Implicit form of the ODE  -------------------- */
101c4762a1bSJed Brown 
102c4762a1bSJed Brown static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
103c4762a1bSJed Brown {
104c4762a1bSJed Brown   User              user = (User)ctx;
105c4762a1bSJed Brown   const PetscScalar *u,*udot;
106c4762a1bSJed Brown   PetscScalar       *f;
107c4762a1bSJed Brown 
108c4762a1bSJed Brown   PetscFunctionBeginUser;
109*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(U,&u));
110*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(Udot,&udot));
111*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(F,&f));
112c4762a1bSJed Brown   f[0] = udot[0] - u[1];
113c4762a1bSJed Brown   f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]);
114*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(U,&u));
115*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(Udot,&udot));
116*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(F,&f));
117c4762a1bSJed Brown   PetscFunctionReturn(0);
118c4762a1bSJed Brown }
119c4762a1bSJed Brown 
120c4762a1bSJed Brown static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
121c4762a1bSJed Brown {
122c4762a1bSJed Brown   User              user     = (User)ctx;
123c4762a1bSJed Brown   PetscInt          rowcol[] = {0,1};
124c4762a1bSJed Brown   PetscScalar       J[2][2];
125c4762a1bSJed Brown   const PetscScalar *u;
126c4762a1bSJed Brown 
127c4762a1bSJed Brown   PetscFunctionBeginUser;
128*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(U,&u));
129c4762a1bSJed Brown 
130c4762a1bSJed Brown   J[0][0] = a;     J[0][1] =  -1.0;
131c4762a1bSJed Brown   J[1][0] = user->mu*(2.0*u[0]*u[1] + 1.0);   J[1][1] = a - user->mu*(1.0-u[0]*u[0]);
132c4762a1bSJed Brown 
133*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES));
134*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(U,&u));
135c4762a1bSJed Brown 
136*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
137*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
138c4762a1bSJed Brown   if (B && A != B) {
139*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
140*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
141c4762a1bSJed Brown   }
142c4762a1bSJed Brown   PetscFunctionReturn(0);
143c4762a1bSJed Brown }
144c4762a1bSJed Brown 
145c4762a1bSJed Brown static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
146c4762a1bSJed Brown {
147c4762a1bSJed Brown   PetscInt          row[] = {0,1},col[]={0};
148c4762a1bSJed Brown   PetscScalar       J[2][1];
149c4762a1bSJed Brown   const PetscScalar *u;
150c4762a1bSJed Brown 
151c4762a1bSJed Brown   PetscFunctionBeginUser;
152*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(U,&u));
153c4762a1bSJed Brown   J[0][0] = 0;
154c4762a1bSJed Brown   J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
155*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES));
156*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
157*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
158*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(U,&u));
159c4762a1bSJed Brown   PetscFunctionReturn(0);
160c4762a1bSJed Brown }
161c4762a1bSJed Brown 
162c4762a1bSJed Brown /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
163c4762a1bSJed Brown static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
164c4762a1bSJed Brown {
165c4762a1bSJed Brown   PetscErrorCode    ierr;
166c4762a1bSJed Brown   const PetscScalar *u;
167c4762a1bSJed Brown   PetscReal         tfinal, dt;
168c4762a1bSJed Brown   User              user = (User)ctx;
169c4762a1bSJed Brown   Vec               interpolatedU;
170c4762a1bSJed Brown 
171c4762a1bSJed Brown   PetscFunctionBeginUser;
172*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSGetTimeStep(ts,&dt));
173*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSGetMaxTime(ts,&tfinal));
174c4762a1bSJed Brown 
175c4762a1bSJed Brown   while (user->next_output <= t && user->next_output <= tfinal) {
176*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDuplicate(U,&interpolatedU));
177*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSInterpolate(ts,user->next_output,interpolatedU));
178*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(interpolatedU,&u));
179c4762a1bSJed Brown     ierr = PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
180c4762a1bSJed Brown                        (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]),
181c4762a1bSJed Brown                        (double)PetscRealPart(u[1]));CHKERRQ(ierr);
182*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(interpolatedU,&u));
183*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&interpolatedU));
184c4762a1bSJed Brown     user->next_output += 0.1;
185c4762a1bSJed Brown   }
186c4762a1bSJed Brown   PetscFunctionReturn(0);
187c4762a1bSJed Brown }
188c4762a1bSJed Brown 
189c4762a1bSJed Brown int main(int argc,char **argv)
190c4762a1bSJed Brown {
191c4762a1bSJed Brown   TS             ts;
192c4762a1bSJed Brown   PetscBool      monitor = PETSC_FALSE,implicitform = PETSC_TRUE;
193c4762a1bSJed Brown   PetscScalar    *x_ptr,*y_ptr,derp;
194c4762a1bSJed Brown   PetscMPIInt    size;
195c4762a1bSJed Brown   struct _n_User user;
196c4762a1bSJed Brown   PetscErrorCode ierr;
197c4762a1bSJed Brown 
198c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199c4762a1bSJed Brown      Initialize program
200c4762a1bSJed Brown      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201c4762a1bSJed Brown   ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
202*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
2033c633725SBarry Smith   PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
204c4762a1bSJed Brown 
205c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206c4762a1bSJed Brown     Set runtime options
207c4762a1bSJed Brown     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208c4762a1bSJed Brown   user.next_output = 0.0;
209c4762a1bSJed Brown   user.mu          = 1.0e3;
210c4762a1bSJed Brown   user.steps       = 0;
211c4762a1bSJed Brown   user.ftime       = 0.5;
212*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL));
213*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL));
214*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL));
215c4762a1bSJed Brown 
216c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217c4762a1bSJed Brown     Create necessary matrix and vectors, solve same ODE on every process
218c4762a1bSJed Brown     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
219*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.A));
220*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2));
221*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetFromOptions(user.A));
222*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(user.A));
223*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreateVecs(user.A,&user.U,NULL));
224c4762a1bSJed Brown 
225*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jacp));
226*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1));
227*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetFromOptions(user.Jacp));
228*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(user.Jacp));
229c4762a1bSJed Brown 
230c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231c4762a1bSJed Brown      Create timestepping solver context
232c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts));
234*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
235c4762a1bSJed Brown   if (implicitform) {
236*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSSetIFunction(ts,NULL,IFunction,&user));
237*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSSetIJacobian(ts,user.A,user.A,IJacobian,&user));
238*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSSetType(ts,TSCN));
239c4762a1bSJed Brown   } else {
240*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,&user));
241*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSSetRHSJacobian(ts,user.A,user.A,RHSJacobian,&user));
242*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSSetType(ts,TSRK));
243c4762a1bSJed Brown   }
244*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user));
245*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetMaxTime(ts,user.ftime));
246*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetTimeStep(ts,0.001));
247*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP));
248c4762a1bSJed Brown   if (monitor) {
249*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TSMonitorSet(ts,Monitor,&user,NULL));
250c4762a1bSJed Brown   }
251c4762a1bSJed Brown 
252c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253c4762a1bSJed Brown      Set initial conditions
254c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.U,&x_ptr));
256c4762a1bSJed Brown   x_ptr[0] = 2.0;
257c4762a1bSJed Brown   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
258*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.U,&x_ptr));
259*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetTimeStep(ts,0.001));
260c4762a1bSJed Brown 
261c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
262c4762a1bSJed Brown     Save trajectory of solution so that TSAdjointSolve() may be used
263c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
264*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetSaveTrajectory(ts));
265c4762a1bSJed Brown 
266c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267c4762a1bSJed Brown      Set runtime options
268c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
269*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetFromOptions(ts));
270c4762a1bSJed Brown 
271*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSolve(ts,user.U));
272*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSGetSolveTime(ts,&user.ftime));
273*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSGetStepNumber(ts,&user.steps));
274c4762a1bSJed Brown 
275c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276c4762a1bSJed Brown      Adjoint model starts here
277c4762a1bSJed Brown      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
278*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreateVecs(user.A,&user.lambda[0],NULL));
279c4762a1bSJed Brown   /* Set initial conditions for the adjoint integration */
280*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.lambda[0],&y_ptr));
281c4762a1bSJed Brown   y_ptr[0] = 1.0; y_ptr[1] = 0.0;
282*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.lambda[0],&y_ptr));
283*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreateVecs(user.A,&user.lambda[1],NULL));
284*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.lambda[1],&y_ptr));
285c4762a1bSJed Brown   y_ptr[0] = 0.0; y_ptr[1] = 1.0;
286*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.lambda[1],&y_ptr));
287c4762a1bSJed Brown 
288*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreateVecs(user.Jacp,&user.mup[0],NULL));
289*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.mup[0],&x_ptr));
290c4762a1bSJed Brown   x_ptr[0] = 0.0;
291*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.mup[0],&x_ptr));
292*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreateVecs(user.Jacp,&user.mup[1],NULL));
293*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.mup[1],&x_ptr));
294c4762a1bSJed Brown   x_ptr[0] = 0.0;
295*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.mup[1],&x_ptr));
296c4762a1bSJed Brown 
297*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetCostGradients(ts,2,user.lambda,user.mup));
298c4762a1bSJed Brown 
299*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSAdjointSolve(ts));
300c4762a1bSJed Brown 
301*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[y(tf)]/d[y0]  d[y(tf)]/d[z0]\n"));
302*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecView(user.lambda[0],PETSC_VIEWER_STDOUT_WORLD));
303*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[z(tf)]/d[y0]  d[z(tf)]/d[z0]\n"));
304*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecView(user.lambda[1],PETSC_VIEWER_STDOUT_WORLD));
305c4762a1bSJed Brown 
306*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.mup[0],&x_ptr));
307*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.lambda[0],&y_ptr));
308c4762a1bSJed Brown   derp = y_ptr[1]*(-10.0/(81.0*user.mu*user.mu)+2.0*292.0/(2187.0*user.mu*user.mu*user.mu))+x_ptr[0];
309*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.mup[0],&x_ptr));
310*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.lambda[0],&y_ptr));
311*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameters: d[y(tf)]/d[mu]\n%g\n",(double)PetscRealPart(derp)));
312c4762a1bSJed Brown 
313*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.mup[1],&x_ptr));
314*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(user.lambda[1],&y_ptr));
315c4762a1bSJed Brown   derp = y_ptr[1]*(-10.0/(81.0*user.mu*user.mu)+2.0*292.0/(2187.0*user.mu*user.mu*user.mu))+x_ptr[0];
316*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.mup[1],&x_ptr));
317*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(user.lambda[1],&y_ptr));
318*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n sensivitity wrt parameters: d[z(tf)]/d[mu]\n%g\n",(double)PetscRealPart(derp)));
319c4762a1bSJed Brown 
320c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
321c4762a1bSJed Brown      Free work space.  All PETSc objects should be destroyed when they
322c4762a1bSJed Brown      are no longer needed.
323c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
324*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&user.A));
325*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&user.Jacp));
326*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.U));
327*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.lambda[0]));
328*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.lambda[1]));
329*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.mup[0]));
330*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.mup[1]));
331*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSDestroy(&ts));
332c4762a1bSJed Brown 
333c4762a1bSJed Brown   ierr = PetscFinalize();
334c4762a1bSJed Brown   return(ierr);
335c4762a1bSJed Brown }
336c4762a1bSJed Brown 
337c4762a1bSJed Brown /*TEST
338c4762a1bSJed Brown 
339c4762a1bSJed Brown     test:
340c4762a1bSJed Brown       requires: revolve
341c4762a1bSJed Brown       args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -viewer_binary_skip_info -ts_dt 0.001 -mu 100000
342c4762a1bSJed Brown 
343c4762a1bSJed Brown     test:
344c4762a1bSJed Brown       suffix: 2
345c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only
346c4762a1bSJed Brown 
347c4762a1bSJed Brown     test:
348c4762a1bSJed Brown       suffix: 3
349c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only 0
350c4762a1bSJed Brown       output_file: output/ex20adj_2.out
351c4762a1bSJed Brown 
352c4762a1bSJed Brown     test:
353c4762a1bSJed Brown       suffix: 4
354c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack
355c4762a1bSJed Brown       output_file: output/ex20adj_2.out
356c4762a1bSJed Brown 
357c4762a1bSJed Brown     test:
358c4762a1bSJed Brown       suffix: 5
359c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack
360c4762a1bSJed Brown       output_file: output/ex20adj_2.out
361c4762a1bSJed Brown 
362c4762a1bSJed Brown     test:
363c4762a1bSJed Brown       suffix: 6
364c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack 0
365c4762a1bSJed Brown       output_file: output/ex20adj_2.out
366c4762a1bSJed Brown 
367c4762a1bSJed Brown     test:
368c4762a1bSJed Brown       suffix: 7
369c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack 0
370c4762a1bSJed Brown       output_file: output/ex20adj_2.out
371c4762a1bSJed Brown 
372c4762a1bSJed Brown     test:
373c4762a1bSJed Brown       suffix: 8
37460f0b76eSHong Zhang       requires: revolve !cams
375c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_solution_only -ts_trajectory_monitor
376c4762a1bSJed Brown       output_file: output/ex20adj_3.out
377c4762a1bSJed Brown 
378c4762a1bSJed Brown     test:
379c4762a1bSJed Brown       suffix: 9
38060f0b76eSHong Zhang       requires: revolve !cams
381c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_solution_only 0 -ts_trajectory_monitor
382c4762a1bSJed Brown       output_file: output/ex20adj_4.out
383c4762a1bSJed Brown 
384c4762a1bSJed Brown     test:
385c4762a1bSJed Brown       requires: revolve
386c4762a1bSJed Brown       suffix: 10
387c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_revolve_online -ts_trajectory_solution_only
388c4762a1bSJed Brown       output_file: output/ex20adj_2.out
389c4762a1bSJed Brown 
390c4762a1bSJed Brown     test:
391c4762a1bSJed Brown       requires: revolve
392c4762a1bSJed Brown       suffix: 11
393c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_revolve_online -ts_trajectory_solution_only 0
394c4762a1bSJed Brown       output_file: output/ex20adj_2.out
395c4762a1bSJed Brown 
396c4762a1bSJed Brown     test:
397c4762a1bSJed Brown       suffix: 12
398c4762a1bSJed Brown       requires: revolve
399c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_solution_only
400c4762a1bSJed Brown       output_file: output/ex20adj_2.out
401c4762a1bSJed Brown 
402c4762a1bSJed Brown     test:
403c4762a1bSJed Brown       suffix: 13
404c4762a1bSJed Brown       requires: revolve
405c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_solution_only 0
406c4762a1bSJed Brown       output_file: output/ex20adj_2.out
407c4762a1bSJed Brown 
408c4762a1bSJed Brown     test:
409c4762a1bSJed Brown       suffix: 14
410c4762a1bSJed Brown       requires: revolve
411c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack
412c4762a1bSJed Brown       output_file: output/ex20adj_2.out
413c4762a1bSJed Brown 
414c4762a1bSJed Brown     test:
415c4762a1bSJed Brown       suffix: 15
416c4762a1bSJed Brown       requires: revolve
417c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack 0
418c4762a1bSJed Brown       output_file: output/ex20adj_2.out
419c4762a1bSJed Brown 
420c4762a1bSJed Brown     test:
421c4762a1bSJed Brown       suffix: 16
422c4762a1bSJed Brown       requires: revolve
423c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack
424c4762a1bSJed Brown       output_file: output/ex20adj_2.out
425c4762a1bSJed Brown 
426c4762a1bSJed Brown     test:
427c4762a1bSJed Brown       suffix: 17
428c4762a1bSJed Brown       requires: revolve
429c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack 0
430c4762a1bSJed Brown       output_file: output/ex20adj_2.out
431c4762a1bSJed Brown 
432c4762a1bSJed Brown     test:
433c4762a1bSJed Brown       suffix: 18
434c4762a1bSJed Brown       requires: revolve
435c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack
436c4762a1bSJed Brown       output_file: output/ex20adj_2.out
437c4762a1bSJed Brown 
438c4762a1bSJed Brown     test:
439c4762a1bSJed Brown       suffix: 19
440c4762a1bSJed Brown       requires: revolve
441c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack
442c4762a1bSJed Brown       output_file: output/ex20adj_2.out
443c4762a1bSJed Brown 
444c4762a1bSJed Brown     test:
445c4762a1bSJed Brown       suffix: 20
446c4762a1bSJed Brown       requires: revolve
447c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_solution_only 0
448c4762a1bSJed Brown       output_file: output/ex20adj_2.out
449c4762a1bSJed Brown 
450c4762a1bSJed Brown     test:
451c4762a1bSJed Brown       suffix: 21
452c4762a1bSJed Brown       requires: revolve
453c4762a1bSJed Brown       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack 0
454c4762a1bSJed Brown       output_file: output/ex20adj_2.out
455c4762a1bSJed Brown 
456c4762a1bSJed Brown     test:
457c4762a1bSJed Brown       suffix: 22
458c4762a1bSJed Brown       args: -ts_type beuler -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only
459c4762a1bSJed Brown       output_file: output/ex20adj_2.out
460cc4f23bcSHong Zhang 
461cc4f23bcSHong Zhang     test:
462cc4f23bcSHong Zhang       suffix: 23
463cc4f23bcSHong Zhang       requires: cams
464533251d3SHong Zhang       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_units_ram 5 -ts_trajectory_solution_only -ts_trajectory_monitor -ts_trajectory_memory_type cams
46560f0b76eSHong Zhang       output_file: output/ex20adj_5.out
466cc4f23bcSHong Zhang 
467cc4f23bcSHong Zhang     test:
468cc4f23bcSHong Zhang       suffix: 24
469cc4f23bcSHong Zhang       requires: cams
470533251d3SHong Zhang       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_units_ram 5 -ts_trajectory_solution_only 0 -ts_trajectory_monitor -ts_trajectory_memory_type cams
47160f0b76eSHong Zhang       output_file: output/ex20adj_6.out
472cc4f23bcSHong Zhang 
473c4762a1bSJed Brown TEST*/
474