1 #include <petsctao.h> 2 #include <petscts.h> 3 4 typedef struct _n_aircraft *Aircraft; 5 struct _n_aircraft { 6 TS ts,quadts; 7 Vec V,W; /* control variables V and W */ 8 PetscInt nsteps; /* number of time steps */ 9 PetscReal ftime; 10 Mat A,H; 11 Mat Jacp,DRDU,DRDP; 12 Vec U,Lambda[1],Mup[1],Lambda2[1],Mup2[1],Dir; 13 Vec rhshp1[1],rhshp2[1],rhshp3[1],rhshp4[1],inthp1[1],inthp2[1],inthp3[1],inthp4[1]; 14 PetscReal lv,lw; 15 PetscBool mf,eh; 16 }; 17 18 PetscErrorCode FormObjFunctionGradient(Tao,Vec,PetscReal *,Vec,void *); 19 PetscErrorCode FormObjHessian(Tao,Vec,Mat,Mat,void *); 20 PetscErrorCode ComputeObjHessianWithSOA(Vec,PetscScalar[],Aircraft); 21 PetscErrorCode MatrixFreeObjHessian(Tao,Vec,Mat,Mat,void *); 22 PetscErrorCode MyMatMult(Mat,Vec,Vec); 23 24 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 25 { 26 Aircraft actx = (Aircraft)ctx; 27 const PetscScalar *u,*v,*w; 28 PetscScalar *f; 29 PetscInt step; 30 PetscErrorCode ierr; 31 32 PetscFunctionBeginUser; 33 ierr = TSGetStepNumber(ts,&step);CHKERRQ(ierr); 34 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 35 ierr = VecGetArrayRead(actx->V,&v);CHKERRQ(ierr); 36 ierr = VecGetArrayRead(actx->W,&w);CHKERRQ(ierr); 37 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 38 f[0] = v[step]*PetscCosReal(w[step]); 39 f[1] = v[step]*PetscSinReal(w[step]); 40 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 41 ierr = VecRestoreArrayRead(actx->V,&v);CHKERRQ(ierr); 42 ierr = VecRestoreArrayRead(actx->W,&w);CHKERRQ(ierr); 43 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx) 48 { 49 Aircraft actx = (Aircraft)ctx; 50 const PetscScalar *u,*v,*w; 51 PetscInt step,rows[2] = {0,1},rowcol[2]; 52 PetscScalar Jp[2][2]; 53 PetscErrorCode ierr; 54 55 PetscFunctionBeginUser; 56 ierr = MatZeroEntries(A);CHKERRQ(ierr); 57 ierr = TSGetStepNumber(ts,&step);CHKERRQ(ierr); 58 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 59 ierr = VecGetArrayRead(actx->V,&v);CHKERRQ(ierr); 60 ierr = VecGetArrayRead(actx->W,&w);CHKERRQ(ierr); 61 62 Jp[0][0] = PetscCosReal(w[step-1]); 63 Jp[0][1] = -v[step-1]*PetscSinReal(w[step-1]); 64 Jp[1][0] = PetscSinReal(w[step-1]); 65 Jp[1][1] = v[step-1]*PetscCosReal(w[step-1]); 66 67 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 68 ierr = VecRestoreArrayRead(actx->V,&v);CHKERRQ(ierr); 69 ierr = VecRestoreArrayRead(actx->W,&w);CHKERRQ(ierr); 70 71 rowcol[0] = 2*(step-1); 72 rowcol[1] = 2*(step-1)+1; 73 ierr = MatSetValues(A,2,rows,2,rowcol,&Jp[0][0],INSERT_VALUES);CHKERRQ(ierr); 74 75 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 76 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 81 { 82 PetscFunctionBeginUser; 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 87 { 88 PetscFunctionBeginUser; 89 PetscFunctionReturn(0); 90 } 91 92 static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 93 { 94 PetscFunctionBeginUser; 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 99 { 100 Aircraft actx = (Aircraft)ctx; 101 const PetscScalar *v,*w,*vl,*vr,*u; 102 PetscScalar *vhv; 103 PetscScalar dJpdP[2][2][2]={{{0}}}; 104 PetscInt step,i,j,k; 105 PetscErrorCode ierr; 106 107 PetscFunctionBeginUser; 108 ierr = TSGetStepNumber(ts,&step);CHKERRQ(ierr); 109 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 110 ierr = VecGetArrayRead(actx->V,&v);CHKERRQ(ierr); 111 ierr = VecGetArrayRead(actx->W,&w);CHKERRQ(ierr); 112 ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr); 113 ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr); 114 ierr = VecSet(VHV[0],0.0);CHKERRQ(ierr); 115 ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr); 116 117 dJpdP[0][0][1] = -PetscSinReal(w[step-1]); 118 dJpdP[0][1][0] = -PetscSinReal(w[step-1]); 119 dJpdP[0][1][1] = -v[step-1]*PetscCosReal(w[step-1]); 120 dJpdP[1][0][1] = PetscCosReal(w[step-1]); 121 dJpdP[1][1][0] = PetscCosReal(w[step-1]); 122 dJpdP[1][1][1] = -v[step-1]*PetscSinReal(w[step-1]); 123 124 for (j=0; j<2; j++) { 125 vhv[2*(step-1)+j] = 0; 126 for (k=0; k<2; k++) 127 for (i=0; i<2; i++) 128 vhv[2*(step-1)+j] += vl[i]*dJpdP[i][j][k]*vr[2*(step-1)+k]; 129 } 130 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 131 ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr); 132 ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr); 133 ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 /* Vl in NULL,updates to VHV must be added */ 138 static PetscErrorCode IntegrandHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 139 { 140 Aircraft actx = (Aircraft)ctx; 141 const PetscScalar *v,*w,*vr,*u; 142 PetscScalar *vhv; 143 PetscScalar dRudU[2][2]={{0}}; 144 PetscInt step,j,k; 145 PetscErrorCode ierr; 146 147 PetscFunctionBeginUser; 148 ierr = TSGetStepNumber(ts,&step);CHKERRQ(ierr); 149 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 150 ierr = VecGetArrayRead(actx->V,&v);CHKERRQ(ierr); 151 ierr = VecGetArrayRead(actx->W,&w);CHKERRQ(ierr); 152 ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr); 153 ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr); 154 155 dRudU[0][0] = 2.0; 156 dRudU[1][1] = 2.0; 157 158 for (j=0; j<2; j++) { 159 vhv[j] = 0; 160 for (k=0; k<2; k++) 161 vhv[j] += dRudU[j][k]*vr[k]; 162 } 163 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 164 ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr); 165 ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 169 static PetscErrorCode IntegrandHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 170 { 171 PetscFunctionBeginUser; 172 PetscFunctionReturn(0); 173 } 174 175 static PetscErrorCode IntegrandHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 176 { 177 PetscFunctionBeginUser; 178 PetscFunctionReturn(0); 179 } 180 181 static PetscErrorCode IntegrandHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 182 { 183 PetscFunctionBeginUser; 184 PetscFunctionReturn(0); 185 } 186 187 188 static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,void *ctx) 189 { 190 Aircraft actx = (Aircraft)ctx; 191 PetscScalar *r; 192 PetscReal dx,dy; 193 const PetscScalar *u; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 198 ierr = VecGetArray(R,&r);CHKERRQ(ierr); 199 dx = u[0] - actx->lv*t*PetscCosReal(actx->lw); 200 dy = u[1] - actx->lv*t*PetscSinReal(actx->lw); 201 r[0] = dx*dx+dy*dy; 202 ierr = VecRestoreArray(R,&r);CHKERRQ(ierr); 203 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,void *ctx) 208 { 209 Aircraft actx = (Aircraft)ctx; 210 PetscScalar drdu[2][1]; 211 const PetscScalar *u; 212 PetscReal dx,dy; 213 PetscInt row[] = {0,1},col[] = {0}; 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 218 dx = u[0] - actx->lv*t*PetscCosReal(actx->lw); 219 dy = u[1] - actx->lv*t*PetscSinReal(actx->lw); 220 drdu[0][0] = 2.*dx; 221 drdu[1][0] = 2.*dy; 222 ierr = MatSetValues(DRDU,2,row,1,col,&drdu[0][0],INSERT_VALUES);CHKERRQ(ierr); 223 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 224 ierr = MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 225 ierr = MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,void *ctx) 230 { 231 PetscErrorCode ierr; 232 233 PetscFunctionBegin; 234 ierr = MatZeroEntries(DRDP);CHKERRQ(ierr); 235 ierr = MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 236 ierr = MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 int main(int argc,char **argv) 241 { 242 Vec P,PL,PU; 243 struct _n_aircraft aircraft; 244 PetscMPIInt size; 245 Tao tao; 246 KSP ksp; 247 PC pc; 248 PetscScalar *u,*p; 249 PetscInt i; 250 PetscErrorCode ierr; 251 252 /* Initialize program */ 253 ierr = PetscInitialize(&argc,&argv,NULL,NULL);if(ierr) return ierr; 254 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 255 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 256 257 /* Parameter settings */ 258 aircraft.ftime = 1.; /* time interval in hour */ 259 aircraft.nsteps = 10; /* number of steps */ 260 aircraft.lv = 2.0; /* leader speed in kmph */ 261 aircraft.lw = PETSC_PI/4.; /* leader heading angle */ 262 263 ierr = PetscOptionsGetReal(NULL,NULL,"-ftime",&aircraft.ftime,NULL);CHKERRQ(ierr); 264 ierr = PetscOptionsGetInt(NULL,NULL,"-nsteps",&aircraft.nsteps,NULL);CHKERRQ(ierr); 265 ierr = PetscOptionsHasName(NULL,NULL,"-matrixfree",&aircraft.mf);CHKERRQ(ierr); 266 ierr = PetscOptionsHasName(NULL,NULL,"-exacthessian",&aircraft.eh);CHKERRQ(ierr); 267 268 /* Create TAO solver and set desired solution method */ 269 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 270 ierr = TaoSetType(tao,TAOBQNLS);CHKERRQ(ierr); 271 272 /* Create necessary matrix and vectors, solve same ODE on every process */ 273 ierr = MatCreate(PETSC_COMM_WORLD,&aircraft.A);CHKERRQ(ierr); 274 ierr = MatSetSizes(aircraft.A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 275 ierr = MatSetFromOptions(aircraft.A);CHKERRQ(ierr); 276 ierr = MatSetUp(aircraft.A);CHKERRQ(ierr); 277 ierr = MatAssemblyBegin(aircraft.A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 278 ierr = MatAssemblyEnd(aircraft.A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 279 ierr = MatShift(aircraft.A,1);CHKERRQ(ierr); 280 ierr = MatShift(aircraft.A,-1);CHKERRQ(ierr); 281 282 ierr = MatCreate(PETSC_COMM_WORLD,&aircraft.Jacp);CHKERRQ(ierr); 283 ierr = MatSetSizes(aircraft.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,2*aircraft.nsteps);CHKERRQ(ierr); 284 ierr = MatSetFromOptions(aircraft.Jacp);CHKERRQ(ierr); 285 ierr = MatSetUp(aircraft.Jacp);CHKERRQ(ierr); 286 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,1,NULL,&aircraft.DRDP);CHKERRQ(ierr); 287 ierr = MatSetUp(aircraft.DRDP);CHKERRQ(ierr); 288 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&aircraft.DRDU);CHKERRQ(ierr); 289 ierr = MatSetUp(aircraft.DRDU);CHKERRQ(ierr); 290 291 /* Create timestepping solver context */ 292 ierr = TSCreate(PETSC_COMM_WORLD,&aircraft.ts);CHKERRQ(ierr); 293 ierr = TSSetType(aircraft.ts,TSRK);CHKERRQ(ierr); 294 ierr = TSSetRHSFunction(aircraft.ts,NULL,RHSFunction,&aircraft);CHKERRQ(ierr); 295 ierr = TSSetRHSJacobian(aircraft.ts,aircraft.A,aircraft.A,TSComputeRHSJacobianConstant,&aircraft);CHKERRQ(ierr); 296 ierr = TSSetRHSJacobianP(aircraft.ts,aircraft.Jacp,RHSJacobianP,&aircraft);CHKERRQ(ierr); 297 ierr = TSSetExactFinalTime(aircraft.ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 298 ierr = TSSetEquationType(aircraft.ts,TS_EQ_ODE_EXPLICIT);CHKERRQ(ierr); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 299 300 /* Set initial conditions */ 301 ierr = MatCreateVecs(aircraft.A,&aircraft.U,NULL);CHKERRQ(ierr); 302 ierr = TSSetSolution(aircraft.ts,aircraft.U);CHKERRQ(ierr); 303 ierr = VecGetArray(aircraft.U,&u);CHKERRQ(ierr); 304 u[0] = 1.5; 305 u[1] = 0; 306 ierr = VecRestoreArray(aircraft.U,&u);CHKERRQ(ierr); 307 ierr = VecCreate(PETSC_COMM_WORLD,&aircraft.V);CHKERRQ(ierr); 308 ierr = VecSetSizes(aircraft.V,PETSC_DECIDE,aircraft.nsteps);CHKERRQ(ierr); 309 ierr = VecSetUp(aircraft.V);CHKERRQ(ierr); 310 ierr = VecDuplicate(aircraft.V,&aircraft.W);CHKERRQ(ierr); 311 ierr = VecSet(aircraft.V,1.);CHKERRQ(ierr); 312 ierr = VecSet(aircraft.W,PETSC_PI/4.);CHKERRQ(ierr); 313 314 /* Save trajectory of solution so that TSAdjointSolve() may be used */ 315 ierr = TSSetSaveTrajectory(aircraft.ts);CHKERRQ(ierr); 316 317 /* Set sensitivity context */ 318 ierr = TSCreateQuadratureTS(aircraft.ts,PETSC_FALSE,&aircraft.quadts);CHKERRQ(ierr); 319 ierr = TSSetRHSFunction(aircraft.quadts,NULL,(TSRHSFunction)CostIntegrand,&aircraft);CHKERRQ(ierr); 320 ierr = TSSetRHSJacobian(aircraft.quadts,aircraft.DRDU,aircraft.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&aircraft);CHKERRQ(ierr); 321 ierr = TSSetRHSJacobianP(aircraft.quadts,aircraft.DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&aircraft);CHKERRQ(ierr); 322 ierr = MatCreateVecs(aircraft.A,&aircraft.Lambda[0],NULL);CHKERRQ(ierr); 323 ierr = MatCreateVecs(aircraft.Jacp,&aircraft.Mup[0],NULL);CHKERRQ(ierr); 324 if (aircraft.eh) { 325 ierr = MatCreateVecs(aircraft.A,&aircraft.rhshp1[0],NULL);CHKERRQ(ierr); 326 ierr = MatCreateVecs(aircraft.A,&aircraft.rhshp2[0],NULL);CHKERRQ(ierr); 327 ierr = MatCreateVecs(aircraft.Jacp,&aircraft.rhshp3[0],NULL);CHKERRQ(ierr); 328 ierr = MatCreateVecs(aircraft.Jacp,&aircraft.rhshp4[0],NULL);CHKERRQ(ierr); 329 ierr = MatCreateVecs(aircraft.DRDU,&aircraft.inthp1[0],NULL);CHKERRQ(ierr); 330 ierr = MatCreateVecs(aircraft.DRDU,&aircraft.inthp2[0],NULL);CHKERRQ(ierr); 331 ierr = MatCreateVecs(aircraft.DRDP,&aircraft.inthp3[0],NULL);CHKERRQ(ierr); 332 ierr = MatCreateVecs(aircraft.DRDP,&aircraft.inthp4[0],NULL);CHKERRQ(ierr); 333 ierr = MatCreateVecs(aircraft.Jacp,&aircraft.Dir,NULL);CHKERRQ(ierr); 334 ierr = TSSetRHSHessianProduct(aircraft.ts,aircraft.rhshp1,RHSHessianProductUU,aircraft.rhshp2,RHSHessianProductUP,aircraft.rhshp3,RHSHessianProductPU,aircraft.rhshp4,RHSHessianProductPP,&aircraft);CHKERRQ(ierr); 335 ierr = TSSetRHSHessianProduct(aircraft.quadts,aircraft.inthp1,IntegrandHessianProductUU,aircraft.inthp2,IntegrandHessianProductUP,aircraft.inthp3,IntegrandHessianProductPU,aircraft.inthp4,IntegrandHessianProductPP,&aircraft);CHKERRQ(ierr); 336 ierr = MatCreateVecs(aircraft.A,&aircraft.Lambda2[0],NULL);CHKERRQ(ierr); 337 ierr = MatCreateVecs(aircraft.Jacp,&aircraft.Mup2[0],NULL);CHKERRQ(ierr); 338 } 339 ierr = TSSetFromOptions(aircraft.ts);CHKERRQ(ierr); 340 ierr = TSSetMaxTime(aircraft.ts,aircraft.ftime);CHKERRQ(ierr); 341 ierr = TSSetTimeStep(aircraft.ts,aircraft.ftime/aircraft.nsteps);CHKERRQ(ierr); 342 343 /* Set initial solution guess */ 344 ierr = MatCreateVecs(aircraft.Jacp,&P,NULL);CHKERRQ(ierr); 345 ierr = VecGetArray(P,&p);CHKERRQ(ierr); 346 for (i=0; i<aircraft.nsteps; i++) { 347 p[2*i] = 2.0; 348 p[2*i+1] = PETSC_PI/2.0; 349 } 350 ierr = VecRestoreArray(P,&p);CHKERRQ(ierr); 351 ierr = VecDuplicate(P,&PU);CHKERRQ(ierr); 352 ierr = VecDuplicate(P,&PL);CHKERRQ(ierr); 353 ierr = VecGetArray(PU,&p);CHKERRQ(ierr); 354 for (i=0; i<aircraft.nsteps; i++) { 355 p[2*i] = 2.0; 356 p[2*i+1] = PETSC_PI; 357 } 358 ierr = VecRestoreArray(PU,&p);CHKERRQ(ierr); 359 ierr = VecGetArray(PL,&p);CHKERRQ(ierr); 360 for (i=0; i<aircraft.nsteps; i++) { 361 p[2*i] = 0.0; 362 p[2*i+1] = -PETSC_PI; 363 } 364 ierr = VecRestoreArray(PL,&p);CHKERRQ(ierr); 365 366 ierr = TaoSetInitialVector(tao,P);CHKERRQ(ierr); 367 ierr = TaoSetVariableBounds(tao,PL,PU);CHKERRQ(ierr); 368 /* Set routine for function and gradient evaluation */ 369 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormObjFunctionGradient,(void *)&aircraft);CHKERRQ(ierr); 370 371 if (aircraft.eh) { 372 if (aircraft.mf) { 373 ierr = MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,2*aircraft.nsteps,(void*)&aircraft,&aircraft.H);CHKERRQ(ierr); 374 ierr = MatShellSetOperation(aircraft.H,MATOP_MULT,(void(*)(void))MyMatMult);CHKERRQ(ierr); 375 ierr = MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 376 ierr = TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,MatrixFreeObjHessian,(void*)&aircraft);CHKERRQ(ierr); 377 } else { 378 ierr = MatCreateDense(MPI_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,2*aircraft.nsteps,2*aircraft.nsteps,NULL,&(aircraft.H));CHKERRQ(ierr); 379 ierr = MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 380 ierr = TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,FormObjHessian,(void *)&aircraft);CHKERRQ(ierr); 381 } 382 } 383 384 /* Check for any TAO command line options */ 385 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 386 if (ksp) { 387 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 388 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 389 } 390 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 391 392 ierr = TaoSolve(tao); CHKERRQ(ierr); 393 ierr = VecView(P,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 394 395 /* Free TAO data structures */ 396 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 397 398 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 399 Free work space. All PETSc objects should be destroyed when they 400 are no longer needed. 401 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 402 ierr = TSDestroy(&aircraft.ts);CHKERRQ(ierr); 403 ierr = MatDestroy(&aircraft.A);CHKERRQ(ierr); 404 ierr = VecDestroy(&aircraft.U);CHKERRQ(ierr); 405 ierr = VecDestroy(&aircraft.V);CHKERRQ(ierr); 406 ierr = VecDestroy(&aircraft.W);CHKERRQ(ierr); 407 ierr = VecDestroy(&P);CHKERRQ(ierr); 408 ierr = VecDestroy(&PU);CHKERRQ(ierr); 409 ierr = VecDestroy(&PL);CHKERRQ(ierr); 410 ierr = MatDestroy(&aircraft.Jacp);CHKERRQ(ierr); 411 ierr = MatDestroy(&aircraft.DRDU);CHKERRQ(ierr); 412 ierr = MatDestroy(&aircraft.DRDP);CHKERRQ(ierr); 413 ierr = VecDestroy(&aircraft.Lambda[0]);CHKERRQ(ierr); 414 ierr = VecDestroy(&aircraft.Mup[0]);CHKERRQ(ierr); 415 ierr = VecDestroy(&P);CHKERRQ(ierr); 416 if (aircraft.eh) { 417 ierr = VecDestroy(&aircraft.Lambda2[0]);CHKERRQ(ierr); 418 ierr = VecDestroy(&aircraft.Mup2[0]);CHKERRQ(ierr); 419 ierr = VecDestroy(&aircraft.Dir);CHKERRQ(ierr); 420 ierr = VecDestroy(&aircraft.rhshp1[0]);CHKERRQ(ierr); 421 ierr = VecDestroy(&aircraft.rhshp2[0]);CHKERRQ(ierr); 422 ierr = VecDestroy(&aircraft.rhshp3[0]);CHKERRQ(ierr); 423 ierr = VecDestroy(&aircraft.rhshp4[0]);CHKERRQ(ierr); 424 ierr = VecDestroy(&aircraft.inthp1[0]);CHKERRQ(ierr); 425 ierr = VecDestroy(&aircraft.inthp2[0]);CHKERRQ(ierr); 426 ierr = VecDestroy(&aircraft.inthp3[0]);CHKERRQ(ierr); 427 ierr = VecDestroy(&aircraft.inthp4[0]);CHKERRQ(ierr); 428 ierr = MatDestroy(&aircraft.H);CHKERRQ(ierr); 429 } 430 ierr = PetscFinalize(); 431 return ierr; 432 } 433 434 /* 435 FormObjFunctionGradient - Evaluates the function and corresponding gradient. 436 437 Input Parameters: 438 tao - the Tao context 439 P - the input vector 440 ctx - optional aircraft-defined context, as set by TaoSetObjectiveAndGradientRoutine() 441 442 Output Parameters: 443 f - the newly evaluated function 444 G - the newly evaluated gradient 445 */ 446 PetscErrorCode FormObjFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx) 447 { 448 Aircraft actx = (Aircraft)ctx; 449 TS ts = actx->ts; 450 Vec Q; 451 const PetscScalar *p,*q; 452 PetscScalar *u,*v,*w; 453 PetscInt i; 454 PetscErrorCode ierr; 455 456 PetscFunctionBeginUser; 457 ierr = VecGetArrayRead(P,&p);CHKERRQ(ierr); 458 ierr = VecGetArray(actx->V,&v);CHKERRQ(ierr); 459 ierr = VecGetArray(actx->W,&w);CHKERRQ(ierr); 460 for (i=0; i<actx->nsteps; i++) { 461 v[i] = p[2*i]; 462 w[i] = p[2*i+1]; 463 } 464 ierr = VecRestoreArrayRead(P,&p);CHKERRQ(ierr); 465 ierr = VecRestoreArray(actx->V,&v);CHKERRQ(ierr); 466 ierr = VecRestoreArray(actx->W,&w);CHKERRQ(ierr); 467 468 ierr = TSSetTime(ts,0.0);CHKERRQ(ierr); 469 ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr); 470 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 471 ierr = TSSetTimeStep(ts,actx->ftime/actx->nsteps);CHKERRQ(ierr); 472 473 /* reinitialize system state */ 474 ierr = VecGetArray(actx->U,&u);CHKERRQ(ierr); 475 u[0] = 2.0; 476 u[1] = 0; 477 ierr = VecRestoreArray(actx->U,&u);CHKERRQ(ierr); 478 479 /* reinitialize the integral value */ 480 ierr = TSGetCostIntegral(ts,&Q);CHKERRQ(ierr); 481 ierr = VecSet(Q,0.0);CHKERRQ(ierr); 482 483 ierr = TSSolve(ts,actx->U);CHKERRQ(ierr); 484 485 /* Reset initial conditions for the adjoint integration */ 486 ierr = VecSet(actx->Lambda[0],0.0);CHKERRQ(ierr); 487 ierr = VecSet(actx->Mup[0],0.0);CHKERRQ(ierr); 488 ierr = TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);CHKERRQ(ierr); 489 490 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 491 ierr = VecCopy(actx->Mup[0],G);CHKERRQ(ierr); 492 ierr = TSGetCostIntegral(ts,&Q);CHKERRQ(ierr); 493 ierr = VecGetArrayRead(Q,&q);CHKERRQ(ierr); 494 *f = q[0]; 495 ierr = VecRestoreArrayRead(Q,&q);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 PetscErrorCode FormObjHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx) 500 { 501 Aircraft actx = (Aircraft)ctx; 502 const PetscScalar *p; 503 PetscScalar *harr,*v,*w,one = 1.0; 504 PetscInt ind[1]; 505 PetscInt *cols,i; 506 Vec Dir; 507 PetscErrorCode ierr; 508 509 PetscFunctionBeginUser; 510 /* set up control parameters */ 511 ierr = VecGetArrayRead(P,&p);CHKERRQ(ierr); 512 ierr = VecGetArray(actx->V,&v);CHKERRQ(ierr); 513 ierr = VecGetArray(actx->W,&w);CHKERRQ(ierr); 514 for (i=0; i<actx->nsteps; i++) { 515 v[i] = p[2*i]; 516 w[i] = p[2*i+1]; 517 } 518 ierr = VecRestoreArrayRead(P,&p);CHKERRQ(ierr); 519 ierr = VecRestoreArray(actx->V,&v);CHKERRQ(ierr); 520 ierr = VecRestoreArray(actx->W,&w);CHKERRQ(ierr); 521 522 ierr = PetscMalloc1(2*actx->nsteps,&harr);CHKERRQ(ierr); 523 ierr = PetscMalloc1(2*actx->nsteps,&cols);CHKERRQ(ierr); 524 for (i=0; i<2*actx->nsteps; i++) cols[i] = i; 525 ierr = VecDuplicate(P,&Dir);CHKERRQ(ierr); 526 for (i=0; i<2*actx->nsteps; i++) { 527 ind[0] = i; 528 ierr = VecSet(Dir,0.0);CHKERRQ(ierr); 529 ierr = VecSetValues(Dir,1,ind,&one,INSERT_VALUES);CHKERRQ(ierr); 530 ierr = VecAssemblyBegin(Dir);CHKERRQ(ierr); 531 ierr = VecAssemblyEnd(Dir);CHKERRQ(ierr); 532 ierr = ComputeObjHessianWithSOA(Dir,harr,actx);CHKERRQ(ierr); 533 ierr = MatSetValues(H,1,ind,2*actx->nsteps,cols,harr,INSERT_VALUES);CHKERRQ(ierr); 534 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 535 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 536 if (H != Hpre) { 537 ierr = MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 538 ierr = MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 539 } 540 } 541 ierr = PetscFree(cols);CHKERRQ(ierr); 542 ierr = PetscFree(harr);CHKERRQ(ierr); 543 ierr = VecDestroy(&Dir);CHKERRQ(ierr); 544 PetscFunctionReturn(0); 545 } 546 547 PetscErrorCode MatrixFreeObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx) 548 { 549 Aircraft actx = (Aircraft)ctx; 550 PetscScalar *v,*w; 551 const PetscScalar *p; 552 PetscInt i; 553 PetscErrorCode ierr; 554 555 PetscFunctionBegin; 556 ierr = VecGetArrayRead(P,&p);CHKERRQ(ierr); 557 ierr = VecGetArray(actx->V,&v);CHKERRQ(ierr); 558 ierr = VecGetArray(actx->W,&w);CHKERRQ(ierr); 559 for (i=0; i<actx->nsteps; i++) { 560 v[i] = p[2*i]; 561 w[i] = p[2*i+1]; 562 } 563 ierr = VecRestoreArrayRead(P,&p);CHKERRQ(ierr); 564 ierr = VecRestoreArray(actx->V,&v);CHKERRQ(ierr); 565 ierr = VecRestoreArray(actx->W,&w);CHKERRQ(ierr); 566 PetscFunctionReturn(0); 567 } 568 569 PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y) 570 { 571 PetscScalar *y; 572 void *ptr; 573 PetscErrorCode ierr; 574 575 PetscFunctionBegin; 576 ierr = MatShellGetContext(H_shell,&ptr);CHKERRQ(ierr); 577 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 578 ierr = ComputeObjHessianWithSOA(X,y,(Aircraft)ptr);CHKERRQ(ierr); 579 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 PetscErrorCode ComputeObjHessianWithSOA(Vec Dir,PetscScalar arr[],Aircraft actx) 584 { 585 TS ts = actx->ts; 586 const PetscScalar *z_ptr; 587 PetscScalar *u; 588 Vec Q; 589 PetscInt i; 590 PetscErrorCode ierr; 591 592 PetscFunctionBeginUser; 593 /* Reset TSAdjoint so that AdjointSetUp will be called again */ 594 ierr = TSAdjointReset(ts);CHKERRQ(ierr); 595 596 ierr = TSSetTime(ts,0.0);CHKERRQ(ierr); 597 ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr); 598 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 599 ierr = TSSetTimeStep(ts,actx->ftime/actx->nsteps);CHKERRQ(ierr); 600 ierr = TSSetCostHessianProducts(actx->ts,1,actx->Lambda2,actx->Mup2,Dir);CHKERRQ(ierr); 601 602 /* reinitialize system state */ 603 ierr = VecGetArray(actx->U,&u);CHKERRQ(ierr); 604 u[0] = 2.0; 605 u[1] = 0; 606 ierr = VecRestoreArray(actx->U,&u);CHKERRQ(ierr); 607 608 /* reinitialize the integral value */ 609 ierr = TSGetCostIntegral(ts,&Q);CHKERRQ(ierr); 610 ierr = VecSet(Q,0.0);CHKERRQ(ierr); 611 612 /* initialize tlm variable */ 613 ierr = MatZeroEntries(actx->Jacp);CHKERRQ(ierr); 614 ierr = MatAssemblyBegin(actx->Jacp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 615 ierr = MatAssemblyEnd(actx->Jacp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 616 ierr = TSAdjointSetForward(ts,actx->Jacp);CHKERRQ(ierr); 617 618 ierr = TSSolve(ts,actx->U);CHKERRQ(ierr); 619 620 /* Set terminal conditions for first- and second-order adjonts */ 621 ierr = VecSet(actx->Lambda[0],0.0);CHKERRQ(ierr); 622 ierr = VecSet(actx->Mup[0],0.0);CHKERRQ(ierr); 623 ierr = VecSet(actx->Lambda2[0],0.0);CHKERRQ(ierr); 624 ierr = VecSet(actx->Mup2[0],0.0);CHKERRQ(ierr); 625 ierr = TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);CHKERRQ(ierr); 626 627 ierr = TSGetCostIntegral(ts,&Q);CHKERRQ(ierr); 628 629 /* Reset initial conditions for the adjoint integration */ 630 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 631 632 /* initial condition does not depend on p, so that lambda is not needed to assemble G */ 633 ierr = VecGetArrayRead(actx->Mup2[0],&z_ptr);CHKERRQ(ierr); 634 for (i=0; i<2*actx->nsteps; i++) arr[i] = z_ptr[i]; 635 ierr = VecRestoreArrayRead(actx->Mup2[0],&z_ptr);CHKERRQ(ierr); 636 637 /* Disable second-order adjoint mode */ 638 ierr = TSAdjointReset(ts);CHKERRQ(ierr); 639 ierr = TSAdjointResetForward(ts);CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 /*TEST 644 build: 645 requires: !complex !single 646 647 test: 648 args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_gatol 1e-7 649 650 test: 651 suffix: 2 652 args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_view -tao_type bntr -tao_bnk_pc_type none -exacthessian 653 654 test: 655 suffix: 3 656 args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_view -tao_type bntr -tao_bnk_pc_type none -exacthessian -matrixfree 657 TEST*/ 658 659 660