1 #include <petsctao.h> 2 3 typedef struct { 4 PetscInt n; /* Number of variables */ 5 PetscInt m; /* Number of constraints */ 6 PetscInt mx; /* grid points in each direction */ 7 PetscInt nt; /* Number of time steps */ 8 PetscInt ndata; /* Number of data points per sample */ 9 IS s_is; 10 IS d_is; 11 VecScatter state_scatter; 12 VecScatter design_scatter; 13 VecScatter *uxi_scatter, *uyi_scatter, *ux_scatter, *uy_scatter, *ui_scatter; 14 VecScatter *yi_scatter; 15 16 Mat Js, Jd, JsBlockPrec, JsInv, JsBlock; 17 PetscBool jformed, c_formed; 18 19 PetscReal alpha; /* Regularization parameter */ 20 PetscReal gamma; 21 PetscReal ht; /* Time step */ 22 PetscReal T; /* Final time */ 23 Mat Q, QT; 24 Mat L, LT; 25 Mat Div, Divwork, Divxy[2]; 26 Mat Grad, Gradxy[2]; 27 Mat M; 28 Mat *C, *Cwork; 29 /* Mat Hs,Hd,Hsd; */ 30 Vec q; 31 Vec ur; /* reference */ 32 33 Vec d; 34 Vec dwork; 35 36 Vec y; /* state variables */ 37 Vec ywork; 38 Vec ytrue; 39 Vec *yi, *yiwork, *ziwork; 40 Vec *uxi, *uyi, *uxiwork, *uyiwork, *ui, *uiwork; 41 42 Vec u; /* design variables */ 43 Vec uwork, vwork; 44 Vec utrue; 45 46 Vec js_diag; 47 48 Vec c; /* constraint vector */ 49 Vec cwork; 50 51 Vec lwork; 52 53 KSP solver; 54 PC prec; 55 PetscInt block_index; 56 57 PetscInt ksp_its; 58 PetscInt ksp_its_initial; 59 } AppCtx; 60 61 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *); 62 PetscErrorCode FormGradient(Tao, Vec, Vec, void *); 63 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *); 64 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void *); 65 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void *); 66 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *); 67 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *); 68 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 69 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 70 PetscErrorCode HyperbolicInitialize(AppCtx *user); 71 PetscErrorCode HyperbolicDestroy(AppCtx *user); 72 PetscErrorCode HyperbolicMonitor(Tao, void *); 73 74 PetscErrorCode StateMatMult(Mat, Vec, Vec); 75 PetscErrorCode StateMatBlockMult(Mat, Vec, Vec); 76 PetscErrorCode StateMatBlockMultTranspose(Mat, Vec, Vec); 77 PetscErrorCode StateMatMultTranspose(Mat, Vec, Vec); 78 PetscErrorCode StateMatGetDiagonal(Mat, Vec); 79 PetscErrorCode StateMatDuplicate(Mat, MatDuplicateOption, Mat *); 80 PetscErrorCode StateMatInvMult(Mat, Vec, Vec); 81 PetscErrorCode StateMatInvTransposeMult(Mat, Vec, Vec); 82 PetscErrorCode StateMatBlockPrecMult(PC, Vec, Vec); 83 84 PetscErrorCode DesignMatMult(Mat, Vec, Vec); 85 PetscErrorCode DesignMatMultTranspose(Mat, Vec, Vec); 86 87 PetscErrorCode Scatter_yi(Vec, Vec *, VecScatter *, PetscInt); /* y to y1,y2,...,y_nt */ 88 PetscErrorCode Gather_yi(Vec, Vec *, VecScatter *, PetscInt); 89 PetscErrorCode Scatter_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt); /* u to ux_1,uy_1,ux_2,uy_2,...,u */ 90 PetscErrorCode Gather_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt); 91 92 static char help[] = ""; 93 94 int main(int argc, char **argv) 95 { 96 Vec x, x0; 97 Tao tao; 98 AppCtx user; 99 IS is_allstate, is_alldesign; 100 PetscInt lo, hi, hi2, lo2, ksp_old; 101 PetscInt ntests = 1; 102 PetscInt i; 103 PetscLogStage stages[1]; 104 105 PetscFunctionBeginUser; 106 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 107 user.mx = 32; 108 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "hyperbolic example", NULL); 109 PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL)); 110 user.nt = 16; 111 PetscCall(PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL)); 112 user.ndata = 64; 113 PetscCall(PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.ndata, NULL)); 114 user.alpha = 10.0; 115 PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL)); 116 user.T = 1.0 / 32.0; 117 PetscCall(PetscOptionsReal("-Tfinal", "Final time", "", user.T, &user.T, NULL)); 118 PetscCall(PetscOptionsInt("-ntests", "Number of times to repeat TaoSolve", "", ntests, &ntests, NULL)); 119 PetscOptionsEnd(); 120 121 user.m = user.mx * user.mx * user.nt; /* number of constraints */ 122 user.n = user.mx * user.mx * 3 * user.nt; /* number of variables */ 123 user.ht = user.T / user.nt; /* Time step */ 124 user.gamma = user.T * user.ht / (user.mx * user.mx); 125 126 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u)); 127 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y)); 128 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c)); 129 PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m)); 130 PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m)); 131 PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m)); 132 PetscCall(VecSetFromOptions(user.u)); 133 PetscCall(VecSetFromOptions(user.y)); 134 PetscCall(VecSetFromOptions(user.c)); 135 136 /* Create scatters for reduced spaces. 137 If the state vector y and design vector u are partitioned as 138 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors), 139 then the solution vector x is organized as 140 [y_1; u_1; y_2; u_2; ...; y_np; u_np]. 141 The index sets user.s_is and user.d_is correspond to the indices of the 142 state and design variables owned by the current processor. 143 */ 144 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 145 146 PetscCall(VecGetOwnershipRange(user.y, &lo, &hi)); 147 PetscCall(VecGetOwnershipRange(user.u, &lo2, &hi2)); 148 149 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_allstate)); 150 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + lo2, 1, &user.s_is)); 151 152 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, lo2, 1, &is_alldesign)); 153 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, hi + lo2, 1, &user.d_is)); 154 155 PetscCall(VecSetSizes(x, hi - lo + hi2 - lo2, user.n)); 156 PetscCall(VecSetFromOptions(x)); 157 158 PetscCall(VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter)); 159 PetscCall(VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter)); 160 PetscCall(ISDestroy(&is_alldesign)); 161 PetscCall(ISDestroy(&is_allstate)); 162 163 /* Create TAO solver and set desired solution method */ 164 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 165 PetscCall(TaoSetType(tao, TAOLCL)); 166 167 /* Set up initial vectors and matrices */ 168 PetscCall(HyperbolicInitialize(&user)); 169 170 PetscCall(Gather(x, user.y, user.state_scatter, user.u, user.design_scatter)); 171 PetscCall(VecDuplicate(x, &x0)); 172 PetscCall(VecCopy(x, x0)); 173 174 /* Set solution vector with an initial guess */ 175 PetscCall(TaoSetSolution(tao, x)); 176 PetscCall(TaoSetObjective(tao, FormFunction, &user)); 177 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user)); 178 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user)); 179 PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, &user)); 180 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user)); 181 PetscCall(TaoSetFromOptions(tao)); 182 PetscCall(TaoSetStateDesignIS(tao, user.s_is, user.d_is)); 183 184 /* SOLVE THE APPLICATION */ 185 PetscCall(PetscLogStageRegister("Trials", &stages[0])); 186 PetscCall(PetscLogStagePush(stages[0])); 187 user.ksp_its_initial = user.ksp_its; 188 ksp_old = user.ksp_its; 189 for (i = 0; i < ntests; i++) { 190 PetscCall(TaoSolve(tao)); 191 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp_old)); 192 PetscCall(VecCopy(x0, x)); 193 PetscCall(TaoSetSolution(tao, x)); 194 } 195 PetscCall(PetscLogStagePop()); 196 PetscCall(PetscBarrier((PetscObject)x)); 197 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations within initialization: ")); 198 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial)); 199 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Total KSP iterations over %" PetscInt_FMT " trial(s): ", ntests)); 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its)); 201 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations per trial: ")); 202 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initial) / ntests)); 203 204 PetscCall(TaoDestroy(&tao)); 205 PetscCall(VecDestroy(&x)); 206 PetscCall(VecDestroy(&x0)); 207 PetscCall(HyperbolicDestroy(&user)); 208 PetscCall(PetscFinalize()); 209 return 0; 210 } 211 /* ------------------------------------------------------------------- */ 212 /* 213 dwork = Qy - d 214 lwork = L*(u-ur).^2 215 f = 1/2 * (dwork.dork + alpha*y.lwork) 216 */ 217 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr) 218 { 219 PetscReal d1 = 0, d2 = 0; 220 AppCtx *user = (AppCtx *)ptr; 221 222 PetscFunctionBegin; 223 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 224 PetscCall(MatMult(user->Q, user->y, user->dwork)); 225 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 226 PetscCall(VecDot(user->dwork, user->dwork, &d1)); 227 228 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); 229 PetscCall(VecPointwiseMult(user->uwork, user->uwork, user->uwork)); 230 PetscCall(MatMult(user->L, user->uwork, user->lwork)); 231 PetscCall(VecDot(user->y, user->lwork, &d2)); 232 *f = 0.5 * (d1 + user->alpha * d2); 233 PetscFunctionReturn(PETSC_SUCCESS); 234 } 235 236 /* ------------------------------------------------------------------- */ 237 /* 238 state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2 239 design: g_d = alpha*(L'y).*(u-ur) 240 */ 241 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr) 242 { 243 AppCtx *user = (AppCtx *)ptr; 244 245 PetscFunctionBegin; 246 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 247 PetscCall(MatMult(user->Q, user->y, user->dwork)); 248 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 249 250 PetscCall(MatMult(user->QT, user->dwork, user->ywork)); 251 252 PetscCall(MatMult(user->LT, user->y, user->uwork)); 253 PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u)); 254 PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork)); 255 PetscCall(VecScale(user->uwork, user->alpha)); 256 257 PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork)); 258 PetscCall(MatMult(user->L, user->vwork, user->lwork)); 259 PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork)); 260 261 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 262 PetscFunctionReturn(PETSC_SUCCESS); 263 } 264 265 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) 266 { 267 PetscReal d1, d2; 268 AppCtx *user = (AppCtx *)ptr; 269 270 PetscFunctionBegin; 271 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 272 PetscCall(MatMult(user->Q, user->y, user->dwork)); 273 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 274 275 PetscCall(MatMult(user->QT, user->dwork, user->ywork)); 276 277 PetscCall(VecDot(user->dwork, user->dwork, &d1)); 278 279 PetscCall(MatMult(user->LT, user->y, user->uwork)); 280 PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u)); 281 PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork)); 282 PetscCall(VecScale(user->uwork, user->alpha)); 283 284 PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork)); 285 PetscCall(MatMult(user->L, user->vwork, user->lwork)); 286 PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork)); 287 288 PetscCall(VecDot(user->y, user->lwork, &d2)); 289 290 *f = 0.5 * (d1 + user->alpha * d2); 291 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 295 /* ------------------------------------------------------------------- */ 296 /* A 297 MatShell object 298 */ 299 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr) 300 { 301 PetscInt i; 302 AppCtx *user = (AppCtx *)ptr; 303 304 PetscFunctionBegin; 305 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 306 PetscCall(Scatter_yi(user->u, user->ui, user->ui_scatter, user->nt)); 307 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 308 for (i = 0; i < user->nt; i++) { 309 PetscCall(MatCopy(user->Divxy[0], user->C[i], SUBSET_NONZERO_PATTERN)); 310 PetscCall(MatCopy(user->Divxy[1], user->Cwork[i], SAME_NONZERO_PATTERN)); 311 312 PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i])); 313 PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i])); 314 PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], SUBSET_NONZERO_PATTERN)); 315 PetscCall(MatScale(user->C[i], user->ht)); 316 PetscCall(MatShift(user->C[i], 1.0)); 317 } 318 PetscFunctionReturn(PETSC_SUCCESS); 319 } 320 321 /* ------------------------------------------------------------------- */ 322 /* B */ 323 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr) 324 { 325 AppCtx *user = (AppCtx *)ptr; 326 327 PetscFunctionBegin; 328 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 329 PetscFunctionReturn(PETSC_SUCCESS); 330 } 331 332 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y) 333 { 334 PetscInt i; 335 AppCtx *user; 336 337 PetscFunctionBegin; 338 PetscCall(MatShellGetContext(J_shell, &user)); 339 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 340 user->block_index = 0; 341 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 342 343 for (i = 1; i < user->nt; i++) { 344 user->block_index = i; 345 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 346 PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1])); 347 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1])); 348 } 349 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 350 PetscFunctionReturn(PETSC_SUCCESS); 351 } 352 353 PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y) 354 { 355 PetscInt i; 356 AppCtx *user; 357 358 PetscFunctionBegin; 359 PetscCall(MatShellGetContext(J_shell, &user)); 360 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 361 362 for (i = 0; i < user->nt - 1; i++) { 363 user->block_index = i; 364 PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i])); 365 PetscCall(MatMult(user->M, user->yi[i + 1], user->ziwork[i + 1])); 366 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i + 1])); 367 } 368 369 i = user->nt - 1; 370 user->block_index = i; 371 PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i])); 372 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 373 PetscFunctionReturn(PETSC_SUCCESS); 374 } 375 376 PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y) 377 { 378 PetscInt i; 379 AppCtx *user; 380 381 PetscFunctionBegin; 382 PetscCall(MatShellGetContext(J_shell, &user)); 383 i = user->block_index; 384 PetscCall(VecPointwiseMult(user->uxiwork[i], X, user->uxi[i])); 385 PetscCall(VecPointwiseMult(user->uyiwork[i], X, user->uyi[i])); 386 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 387 PetscCall(MatMult(user->Div, user->uiwork[i], Y)); 388 PetscCall(VecAYPX(Y, user->ht, X)); 389 PetscFunctionReturn(PETSC_SUCCESS); 390 } 391 392 PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y) 393 { 394 PetscInt i; 395 AppCtx *user; 396 397 PetscFunctionBegin; 398 PetscCall(MatShellGetContext(J_shell, &user)); 399 i = user->block_index; 400 PetscCall(MatMult(user->Grad, X, user->uiwork[i])); 401 PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 402 PetscCall(VecPointwiseMult(user->uxiwork[i], user->uxi[i], user->uxiwork[i])); 403 PetscCall(VecPointwiseMult(user->uyiwork[i], user->uyi[i], user->uyiwork[i])); 404 PetscCall(VecWAXPY(Y, 1.0, user->uxiwork[i], user->uyiwork[i])); 405 PetscCall(VecAYPX(Y, user->ht, X)); 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408 409 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y) 410 { 411 PetscInt i; 412 AppCtx *user; 413 414 PetscFunctionBegin; 415 PetscCall(MatShellGetContext(J_shell, &user)); 416 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt)); 417 PetscCall(Scatter_uxi_uyi(X, user->uxiwork, user->uxi_scatter, user->uyiwork, user->uyi_scatter, user->nt)); 418 for (i = 0; i < user->nt; i++) { 419 PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i])); 420 PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i])); 421 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 422 PetscCall(MatMult(user->Div, user->uiwork[i], user->ziwork[i])); 423 PetscCall(VecScale(user->ziwork[i], user->ht)); 424 } 425 PetscCall(Gather_yi(Y, user->ziwork, user->yi_scatter, user->nt)); 426 PetscFunctionReturn(PETSC_SUCCESS); 427 } 428 429 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 430 { 431 PetscInt i; 432 AppCtx *user; 433 434 PetscFunctionBegin; 435 PetscCall(MatShellGetContext(J_shell, &user)); 436 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt)); 437 PetscCall(Scatter_yi(X, user->yiwork, user->yi_scatter, user->nt)); 438 for (i = 0; i < user->nt; i++) { 439 PetscCall(MatMult(user->Grad, user->yiwork[i], user->uiwork[i])); 440 PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 441 PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i])); 442 PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i])); 443 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 444 PetscCall(VecScale(user->uiwork[i], user->ht)); 445 } 446 PetscCall(Gather_yi(Y, user->uiwork, user->ui_scatter, user->nt)); 447 PetscFunctionReturn(PETSC_SUCCESS); 448 } 449 450 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y) 451 { 452 PetscInt i; 453 AppCtx *user; 454 455 PetscFunctionBegin; 456 PetscCall(PCShellGetContext(PC_shell, &user)); 457 i = user->block_index; 458 if (user->c_formed) { 459 PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y)); 460 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed"); 461 PetscFunctionReturn(PETSC_SUCCESS); 462 } 463 464 PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y) 465 { 466 PetscInt i; 467 AppCtx *user; 468 469 PetscFunctionBegin; 470 PetscCall(PCShellGetContext(PC_shell, &user)); 471 472 i = user->block_index; 473 if (user->c_formed) { 474 PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y)); 475 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed"); 476 PetscFunctionReturn(PETSC_SUCCESS); 477 } 478 479 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y) 480 { 481 AppCtx *user; 482 PetscInt its, i; 483 484 PetscFunctionBegin; 485 PetscCall(MatShellGetContext(J_shell, &user)); 486 487 if (Y == user->ytrue) { 488 /* First solve is done using true solution to set up problem */ 489 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, PETSC_DEFAULT, PETSC_DEFAULT)); 490 } else { 491 PetscCall(KSPSetTolerances(user->solver, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 492 } 493 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 494 PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 495 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 496 497 user->block_index = 0; 498 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0])); 499 500 PetscCall(KSPGetIterationNumber(user->solver, &its)); 501 user->ksp_its = user->ksp_its + its; 502 for (i = 1; i < user->nt; i++) { 503 PetscCall(MatMult(user->M, user->yiwork[i - 1], user->ziwork[i - 1])); 504 PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i - 1])); 505 user->block_index = i; 506 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 507 508 PetscCall(KSPGetIterationNumber(user->solver, &its)); 509 user->ksp_its = user->ksp_its + its; 510 } 511 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 512 PetscFunctionReturn(PETSC_SUCCESS); 513 } 514 515 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y) 516 { 517 AppCtx *user; 518 PetscInt its, i; 519 520 PetscFunctionBegin; 521 PetscCall(MatShellGetContext(J_shell, &user)); 522 523 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 524 PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 525 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 526 527 i = user->nt - 1; 528 user->block_index = i; 529 PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i])); 530 531 PetscCall(KSPGetIterationNumber(user->solver, &its)); 532 user->ksp_its = user->ksp_its + its; 533 534 for (i = user->nt - 2; i >= 0; i--) { 535 PetscCall(MatMult(user->M, user->yiwork[i + 1], user->ziwork[i + 1])); 536 PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i + 1])); 537 user->block_index = i; 538 PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i])); 539 540 PetscCall(KSPGetIterationNumber(user->solver, &its)); 541 user->ksp_its = user->ksp_its + its; 542 } 543 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 544 PetscFunctionReturn(PETSC_SUCCESS); 545 } 546 547 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell) 548 { 549 AppCtx *user; 550 551 PetscFunctionBegin; 552 PetscCall(MatShellGetContext(J_shell, &user)); 553 554 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell)); 555 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (void (*)(void))StateMatMult)); 556 PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate)); 557 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose)); 558 PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal)); 559 PetscFunctionReturn(PETSC_SUCCESS); 560 } 561 562 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X) 563 { 564 AppCtx *user; 565 566 PetscFunctionBegin; 567 PetscCall(MatShellGetContext(J_shell, &user)); 568 PetscCall(VecCopy(user->js_diag, X)); 569 PetscFunctionReturn(PETSC_SUCCESS); 570 } 571 572 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr) 573 { 574 /* con = Ay - q, A = [C(u1) 0 0 ... 0; 575 -M C(u2) 0 ... 0; 576 0 -M C(u3) ... 0; 577 ... ; 578 0 ... -M C(u_nt)] 579 C(u) = eye + ht*Div*[diag(u1); diag(u2)] */ 580 PetscInt i; 581 AppCtx *user = (AppCtx *)ptr; 582 583 PetscFunctionBegin; 584 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 585 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt)); 586 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 587 588 user->block_index = 0; 589 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 590 591 for (i = 1; i < user->nt; i++) { 592 user->block_index = i; 593 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 594 PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1])); 595 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1])); 596 } 597 598 PetscCall(Gather_yi(C, user->yiwork, user->yi_scatter, user->nt)); 599 PetscCall(VecAXPY(C, -1.0, user->q)); 600 PetscFunctionReturn(PETSC_SUCCESS); 601 } 602 603 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 604 { 605 PetscFunctionBegin; 606 PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 607 PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 608 PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 609 PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 610 PetscFunctionReturn(PETSC_SUCCESS); 611 } 612 613 PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt) 614 { 615 PetscInt i; 616 617 PetscFunctionBegin; 618 for (i = 0; i < nt; i++) { 619 PetscCall(VecScatterBegin(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD)); 620 PetscCall(VecScatterEnd(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD)); 621 PetscCall(VecScatterBegin(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD)); 622 PetscCall(VecScatterEnd(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD)); 623 } 624 PetscFunctionReturn(PETSC_SUCCESS); 625 } 626 627 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 628 { 629 PetscFunctionBegin; 630 PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 631 PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 632 PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 633 PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 634 PetscFunctionReturn(PETSC_SUCCESS); 635 } 636 637 PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt) 638 { 639 PetscInt i; 640 641 PetscFunctionBegin; 642 for (i = 0; i < nt; i++) { 643 PetscCall(VecScatterBegin(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 644 PetscCall(VecScatterEnd(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 645 PetscCall(VecScatterBegin(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 646 PetscCall(VecScatterEnd(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 647 } 648 PetscFunctionReturn(PETSC_SUCCESS); 649 } 650 651 PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 652 { 653 PetscInt i; 654 655 PetscFunctionBegin; 656 for (i = 0; i < nt; i++) { 657 PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 658 PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 659 } 660 PetscFunctionReturn(PETSC_SUCCESS); 661 } 662 663 PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 664 { 665 PetscInt i; 666 667 PetscFunctionBegin; 668 for (i = 0; i < nt; i++) { 669 PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 670 PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 671 } 672 PetscFunctionReturn(PETSC_SUCCESS); 673 } 674 675 PetscErrorCode HyperbolicInitialize(AppCtx *user) 676 { 677 PetscInt n, i, j, linear_index, istart, iend, iblock, lo, hi; 678 Vec XX, YY, XXwork, YYwork, yi, uxi, ui, bc; 679 PetscReal h, sum; 680 PetscScalar hinv, neg_hinv, quarter = 0.25, one = 1.0, half_hinv, neg_half_hinv; 681 PetscScalar vx, vy, zero = 0.0; 682 IS is_from_y, is_to_yi, is_from_u, is_to_uxi, is_to_uyi; 683 684 PetscFunctionBegin; 685 user->jformed = PETSC_FALSE; 686 user->c_formed = PETSC_FALSE; 687 688 user->ksp_its = 0; 689 user->ksp_its_initial = 0; 690 691 n = user->mx * user->mx; 692 693 h = 1.0 / user->mx; 694 hinv = user->mx; 695 neg_hinv = -hinv; 696 half_hinv = hinv / 2.0; 697 neg_half_hinv = neg_hinv / 2.0; 698 699 /* Generate Grad matrix */ 700 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad)); 701 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, 2 * n, n)); 702 PetscCall(MatSetFromOptions(user->Grad)); 703 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 3, NULL, 3, NULL)); 704 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 3, NULL)); 705 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend)); 706 707 for (i = istart; i < iend; i++) { 708 if (i < n) { 709 iblock = i / user->mx; 710 j = iblock * user->mx + ((i + user->mx - 1) % user->mx); 711 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 712 j = iblock * user->mx + ((i + 1) % user->mx); 713 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 714 } 715 if (i >= n) { 716 j = (i - user->mx) % n; 717 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 718 j = (j + 2 * user->mx) % n; 719 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 720 } 721 } 722 723 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY)); 724 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY)); 725 726 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[0])); 727 PetscCall(MatSetSizes(user->Gradxy[0], PETSC_DECIDE, PETSC_DECIDE, n, n)); 728 PetscCall(MatSetFromOptions(user->Gradxy[0])); 729 PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[0], 3, NULL, 3, NULL)); 730 PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[0], 3, NULL)); 731 PetscCall(MatGetOwnershipRange(user->Gradxy[0], &istart, &iend)); 732 733 for (i = istart; i < iend; i++) { 734 iblock = i / user->mx; 735 j = iblock * user->mx + ((i + user->mx - 1) % user->mx); 736 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 737 j = iblock * user->mx + ((i + 1) % user->mx); 738 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 739 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &i, &zero, INSERT_VALUES)); 740 } 741 PetscCall(MatAssemblyBegin(user->Gradxy[0], MAT_FINAL_ASSEMBLY)); 742 PetscCall(MatAssemblyEnd(user->Gradxy[0], MAT_FINAL_ASSEMBLY)); 743 744 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[1])); 745 PetscCall(MatSetSizes(user->Gradxy[1], PETSC_DECIDE, PETSC_DECIDE, n, n)); 746 PetscCall(MatSetFromOptions(user->Gradxy[1])); 747 PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[1], 3, NULL, 3, NULL)); 748 PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[1], 3, NULL)); 749 PetscCall(MatGetOwnershipRange(user->Gradxy[1], &istart, &iend)); 750 751 for (i = istart; i < iend; i++) { 752 j = (i + n - user->mx) % n; 753 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 754 j = (j + 2 * user->mx) % n; 755 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 756 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &i, &zero, INSERT_VALUES)); 757 } 758 PetscCall(MatAssemblyBegin(user->Gradxy[1], MAT_FINAL_ASSEMBLY)); 759 PetscCall(MatAssemblyEnd(user->Gradxy[1], MAT_FINAL_ASSEMBLY)); 760 761 /* Generate Div matrix */ 762 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div)); 763 PetscCall(MatTranspose(user->Gradxy[0], MAT_INITIAL_MATRIX, &user->Divxy[0])); 764 PetscCall(MatTranspose(user->Gradxy[1], MAT_INITIAL_MATRIX, &user->Divxy[1])); 765 766 /* Off-diagonal averaging matrix */ 767 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->M)); 768 PetscCall(MatSetSizes(user->M, PETSC_DECIDE, PETSC_DECIDE, n, n)); 769 PetscCall(MatSetFromOptions(user->M)); 770 PetscCall(MatMPIAIJSetPreallocation(user->M, 4, NULL, 4, NULL)); 771 PetscCall(MatSeqAIJSetPreallocation(user->M, 4, NULL)); 772 PetscCall(MatGetOwnershipRange(user->M, &istart, &iend)); 773 774 for (i = istart; i < iend; i++) { 775 /* kron(Id,Av) */ 776 iblock = i / user->mx; 777 j = iblock * user->mx + ((i + user->mx - 1) % user->mx); 778 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 779 j = iblock * user->mx + ((i + 1) % user->mx); 780 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 781 782 /* kron(Av,Id) */ 783 j = (i + user->mx) % n; 784 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 785 j = (i + n - user->mx) % n; 786 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 787 } 788 PetscCall(MatAssemblyBegin(user->M, MAT_FINAL_ASSEMBLY)); 789 PetscCall(MatAssemblyEnd(user->M, MAT_FINAL_ASSEMBLY)); 790 791 /* Generate 2D grid */ 792 PetscCall(VecCreate(PETSC_COMM_WORLD, &XX)); 793 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q)); 794 PetscCall(VecSetSizes(XX, PETSC_DECIDE, n)); 795 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt)); 796 PetscCall(VecSetFromOptions(XX)); 797 PetscCall(VecSetFromOptions(user->q)); 798 799 PetscCall(VecDuplicate(XX, &YY)); 800 PetscCall(VecDuplicate(XX, &XXwork)); 801 PetscCall(VecDuplicate(XX, &YYwork)); 802 PetscCall(VecDuplicate(XX, &user->d)); 803 PetscCall(VecDuplicate(XX, &user->dwork)); 804 805 PetscCall(VecGetOwnershipRange(XX, &istart, &iend)); 806 for (linear_index = istart; linear_index < iend; linear_index++) { 807 i = linear_index % user->mx; 808 j = (linear_index - i) / user->mx; 809 vx = h * (i + 0.5); 810 vy = h * (j + 0.5); 811 PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES)); 812 PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES)); 813 } 814 815 PetscCall(VecAssemblyBegin(XX)); 816 PetscCall(VecAssemblyEnd(XX)); 817 PetscCall(VecAssemblyBegin(YY)); 818 PetscCall(VecAssemblyEnd(YY)); 819 820 /* Compute final density function yT 821 yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2)) 822 yT = yT / (h^2*sum(yT)) */ 823 PetscCall(VecCopy(XX, XXwork)); 824 PetscCall(VecCopy(YY, YYwork)); 825 826 PetscCall(VecShift(XXwork, -0.25)); 827 PetscCall(VecShift(YYwork, -0.25)); 828 829 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 830 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 831 832 PetscCall(VecCopy(XXwork, user->dwork)); 833 PetscCall(VecAXPY(user->dwork, 1.0, YYwork)); 834 PetscCall(VecScale(user->dwork, -30.0)); 835 PetscCall(VecExp(user->dwork)); 836 PetscCall(VecCopy(user->dwork, user->d)); 837 838 PetscCall(VecCopy(XX, XXwork)); 839 PetscCall(VecCopy(YY, YYwork)); 840 841 PetscCall(VecShift(XXwork, -0.75)); 842 PetscCall(VecShift(YYwork, -0.75)); 843 844 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 845 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 846 847 PetscCall(VecCopy(XXwork, user->dwork)); 848 PetscCall(VecAXPY(user->dwork, 1.0, YYwork)); 849 PetscCall(VecScale(user->dwork, -30.0)); 850 PetscCall(VecExp(user->dwork)); 851 852 PetscCall(VecAXPY(user->d, 1.0, user->dwork)); 853 PetscCall(VecShift(user->d, 1.0)); 854 PetscCall(VecSum(user->d, &sum)); 855 PetscCall(VecScale(user->d, 1.0 / (h * h * sum))); 856 857 /* Initial conditions of forward problem */ 858 PetscCall(VecDuplicate(XX, &bc)); 859 PetscCall(VecCopy(XX, XXwork)); 860 PetscCall(VecCopy(YY, YYwork)); 861 862 PetscCall(VecShift(XXwork, -0.5)); 863 PetscCall(VecShift(YYwork, -0.5)); 864 865 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 866 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 867 868 PetscCall(VecWAXPY(bc, 1.0, XXwork, YYwork)); 869 PetscCall(VecScale(bc, -50.0)); 870 PetscCall(VecExp(bc)); 871 PetscCall(VecShift(bc, 1.0)); 872 PetscCall(VecSum(bc, &sum)); 873 PetscCall(VecScale(bc, 1.0 / (h * h * sum))); 874 875 /* Create scatter from y to y_1,y_2,...,y_nt */ 876 /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */ 877 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->yi_scatter)); 878 PetscCall(VecCreate(PETSC_COMM_WORLD, &yi)); 879 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx)); 880 PetscCall(VecSetFromOptions(yi)); 881 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi)); 882 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork)); 883 PetscCall(VecDuplicateVecs(yi, user->nt, &user->ziwork)); 884 for (i = 0; i < user->nt; i++) { 885 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi)); 886 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi)); 887 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + i * user->mx * user->mx, 1, &is_from_y)); 888 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i])); 889 PetscCall(ISDestroy(&is_to_yi)); 890 PetscCall(ISDestroy(&is_from_y)); 891 } 892 893 /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */ 894 /* TODO: reorder for better parallelism */ 895 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uxi_scatter)); 896 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uyi_scatter)); 897 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->ux_scatter)); 898 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uy_scatter)); 899 PetscCall(PetscMalloc1(2 * user->nt * user->mx * user->mx, &user->ui_scatter)); 900 PetscCall(VecCreate(PETSC_COMM_WORLD, &uxi)); 901 PetscCall(VecCreate(PETSC_COMM_WORLD, &ui)); 902 PetscCall(VecSetSizes(uxi, PETSC_DECIDE, user->mx * user->mx)); 903 PetscCall(VecSetSizes(ui, PETSC_DECIDE, 2 * user->mx * user->mx)); 904 PetscCall(VecSetFromOptions(uxi)); 905 PetscCall(VecSetFromOptions(ui)); 906 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxi)); 907 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyi)); 908 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxiwork)); 909 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyiwork)); 910 PetscCall(VecDuplicateVecs(ui, user->nt, &user->ui)); 911 PetscCall(VecDuplicateVecs(ui, user->nt, &user->uiwork)); 912 for (i = 0; i < user->nt; i++) { 913 PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi)); 914 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi)); 915 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u)); 916 PetscCall(VecScatterCreate(user->u, is_from_u, user->uxi[i], is_to_uxi, &user->uxi_scatter[i])); 917 918 PetscCall(ISDestroy(&is_to_uxi)); 919 PetscCall(ISDestroy(&is_from_u)); 920 921 PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi)); 922 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi)); 923 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + (2 * i + 1) * user->mx * user->mx, 1, &is_from_u)); 924 PetscCall(VecScatterCreate(user->u, is_from_u, user->uyi[i], is_to_uyi, &user->uyi_scatter[i])); 925 926 PetscCall(ISDestroy(&is_to_uyi)); 927 PetscCall(ISDestroy(&is_from_u)); 928 929 PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi)); 930 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi)); 931 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_from_u)); 932 PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uxi[i], is_to_uxi, &user->ux_scatter[i])); 933 934 PetscCall(ISDestroy(&is_to_uxi)); 935 PetscCall(ISDestroy(&is_from_u)); 936 937 PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi)); 938 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi)); 939 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + user->mx * user->mx, 1, &is_from_u)); 940 PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uyi[i], is_to_uyi, &user->uy_scatter[i])); 941 942 PetscCall(ISDestroy(&is_to_uyi)); 943 PetscCall(ISDestroy(&is_from_u)); 944 945 PetscCall(VecGetOwnershipRange(user->ui[i], &lo, &hi)); 946 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi)); 947 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u)); 948 PetscCall(VecScatterCreate(user->u, is_from_u, user->ui[i], is_to_uxi, &user->ui_scatter[i])); 949 950 PetscCall(ISDestroy(&is_to_uxi)); 951 PetscCall(ISDestroy(&is_from_u)); 952 } 953 954 /* RHS of forward problem */ 955 PetscCall(MatMult(user->M, bc, user->yiwork[0])); 956 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0)); 957 PetscCall(Gather_yi(user->q, user->yiwork, user->yi_scatter, user->nt)); 958 959 /* Compute true velocity field utrue */ 960 PetscCall(VecDuplicate(user->u, &user->utrue)); 961 for (i = 0; i < user->nt; i++) { 962 PetscCall(VecCopy(YY, user->uxi[i])); 963 PetscCall(VecScale(user->uxi[i], 150.0 * i * user->ht)); 964 PetscCall(VecCopy(XX, user->uyi[i])); 965 PetscCall(VecShift(user->uyi[i], -10.0)); 966 PetscCall(VecScale(user->uyi[i], 15.0 * i * user->ht)); 967 } 968 PetscCall(Gather_uxi_uyi(user->utrue, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 969 970 /* Initial guess and reference model */ 971 PetscCall(VecDuplicate(user->utrue, &user->ur)); 972 for (i = 0; i < user->nt; i++) { 973 PetscCall(VecCopy(XX, user->uxi[i])); 974 PetscCall(VecShift(user->uxi[i], i * user->ht)); 975 PetscCall(VecCopy(YY, user->uyi[i])); 976 PetscCall(VecShift(user->uyi[i], -i * user->ht)); 977 } 978 PetscCall(Gather_uxi_uyi(user->ur, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 979 980 /* Generate regularization matrix L */ 981 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->LT)); 982 PetscCall(MatSetSizes(user->LT, PETSC_DECIDE, PETSC_DECIDE, 2 * n * user->nt, n * user->nt)); 983 PetscCall(MatSetFromOptions(user->LT)); 984 PetscCall(MatMPIAIJSetPreallocation(user->LT, 1, NULL, 1, NULL)); 985 PetscCall(MatSeqAIJSetPreallocation(user->LT, 1, NULL)); 986 PetscCall(MatGetOwnershipRange(user->LT, &istart, &iend)); 987 988 for (i = istart; i < iend; i++) { 989 iblock = (i + n) / (2 * n); 990 j = i - iblock * n; 991 PetscCall(MatSetValues(user->LT, 1, &i, 1, &j, &user->gamma, INSERT_VALUES)); 992 } 993 994 PetscCall(MatAssemblyBegin(user->LT, MAT_FINAL_ASSEMBLY)); 995 PetscCall(MatAssemblyEnd(user->LT, MAT_FINAL_ASSEMBLY)); 996 997 PetscCall(MatTranspose(user->LT, MAT_INITIAL_MATRIX, &user->L)); 998 999 /* Build work vectors and matrices */ 1000 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork)); 1001 PetscCall(VecSetType(user->lwork, VECMPI)); 1002 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, user->m)); 1003 PetscCall(VecSetFromOptions(user->lwork)); 1004 1005 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork)); 1006 1007 PetscCall(VecDuplicate(user->y, &user->ywork)); 1008 PetscCall(VecDuplicate(user->u, &user->uwork)); 1009 PetscCall(VecDuplicate(user->u, &user->vwork)); 1010 PetscCall(VecDuplicate(user->u, &user->js_diag)); 1011 PetscCall(VecDuplicate(user->c, &user->cwork)); 1012 1013 /* Create matrix-free shell user->Js for computing A*x */ 1014 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->Js)); 1015 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (void (*)(void))StateMatMult)); 1016 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate)); 1017 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose)); 1018 PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal)); 1019 1020 /* Diagonal blocks of user->Js */ 1021 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlock)); 1022 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (void (*)(void))StateMatBlockMult)); 1023 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockMultTranspose)); 1024 1025 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U, 1026 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular. 1027 This is an SOR preconditioner for user->JsBlock. */ 1028 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlockPrec)); 1029 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (void (*)(void))StateMatBlockPrecMult)); 1030 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockPrecMultTranspose)); 1031 1032 /* Create a matrix-free shell user->Jd for computing B*x */ 1033 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->n - user->m, user, &user->Jd)); 1034 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (void (*)(void))DesignMatMult)); 1035 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (void (*)(void))DesignMatMultTranspose)); 1036 1037 /* User-defined routines for computing user->Js\x and user->Js^T\x*/ 1038 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsInv)); 1039 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (void (*)(void))StateMatInvMult)); 1040 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatInvTransposeMult)); 1041 1042 /* Build matrices for SOR preconditioner */ 1043 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 1044 PetscCall(PetscMalloc1(5 * n, &user->C)); 1045 PetscCall(PetscMalloc1(2 * n, &user->Cwork)); 1046 for (i = 0; i < user->nt; i++) { 1047 PetscCall(MatDuplicate(user->Divxy[0], MAT_COPY_VALUES, &user->C[i])); 1048 PetscCall(MatDuplicate(user->Divxy[1], MAT_COPY_VALUES, &user->Cwork[i])); 1049 1050 PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i])); 1051 PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i])); 1052 PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], DIFFERENT_NONZERO_PATTERN)); 1053 PetscCall(MatScale(user->C[i], user->ht)); 1054 PetscCall(MatShift(user->C[i], 1.0)); 1055 } 1056 1057 /* Solver options and tolerances */ 1058 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver)); 1059 PetscCall(KSPSetType(user->solver, KSPGMRES)); 1060 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec)); 1061 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500)); 1062 /* PetscCall(KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500)); */ 1063 PetscCall(KSPGetPC(user->solver, &user->prec)); 1064 PetscCall(PCSetType(user->prec, PCSHELL)); 1065 1066 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult)); 1067 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMultTranspose)); 1068 PetscCall(PCShellSetContext(user->prec, user)); 1069 1070 /* Compute true state function yt given ut */ 1071 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue)); 1072 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt)); 1073 PetscCall(VecSetFromOptions(user->ytrue)); 1074 user->c_formed = PETSC_TRUE; 1075 PetscCall(VecCopy(user->utrue, user->u)); /* Set u=utrue temporarily for StateMatInv */ 1076 PetscCall(VecSet(user->ytrue, 0.0)); /* Initial guess */ 1077 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue)); 1078 PetscCall(VecCopy(user->ur, user->u)); /* Reset u=ur */ 1079 1080 /* Initial guess y0 for state given u0 */ 1081 PetscCall(StateMatInvMult(user->Js, user->q, user->y)); 1082 1083 /* Data discretization */ 1084 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Q)); 1085 PetscCall(MatSetSizes(user->Q, PETSC_DECIDE, PETSC_DECIDE, user->mx * user->mx, user->m)); 1086 PetscCall(MatSetFromOptions(user->Q)); 1087 PetscCall(MatMPIAIJSetPreallocation(user->Q, 0, NULL, 1, NULL)); 1088 PetscCall(MatSeqAIJSetPreallocation(user->Q, 1, NULL)); 1089 1090 PetscCall(MatGetOwnershipRange(user->Q, &istart, &iend)); 1091 1092 for (i = istart; i < iend; i++) { 1093 j = i + user->m - user->mx * user->mx; 1094 PetscCall(MatSetValues(user->Q, 1, &i, 1, &j, &one, INSERT_VALUES)); 1095 } 1096 1097 PetscCall(MatAssemblyBegin(user->Q, MAT_FINAL_ASSEMBLY)); 1098 PetscCall(MatAssemblyEnd(user->Q, MAT_FINAL_ASSEMBLY)); 1099 1100 PetscCall(MatTranspose(user->Q, MAT_INITIAL_MATRIX, &user->QT)); 1101 1102 PetscCall(VecDestroy(&XX)); 1103 PetscCall(VecDestroy(&YY)); 1104 PetscCall(VecDestroy(&XXwork)); 1105 PetscCall(VecDestroy(&YYwork)); 1106 PetscCall(VecDestroy(&yi)); 1107 PetscCall(VecDestroy(&uxi)); 1108 PetscCall(VecDestroy(&ui)); 1109 PetscCall(VecDestroy(&bc)); 1110 1111 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */ 1112 PetscCall(KSPSetFromOptions(user->solver)); 1113 PetscFunctionReturn(PETSC_SUCCESS); 1114 } 1115 1116 PetscErrorCode HyperbolicDestroy(AppCtx *user) 1117 { 1118 PetscInt i; 1119 1120 PetscFunctionBegin; 1121 PetscCall(MatDestroy(&user->Q)); 1122 PetscCall(MatDestroy(&user->QT)); 1123 PetscCall(MatDestroy(&user->Div)); 1124 PetscCall(MatDestroy(&user->Divwork)); 1125 PetscCall(MatDestroy(&user->Grad)); 1126 PetscCall(MatDestroy(&user->L)); 1127 PetscCall(MatDestroy(&user->LT)); 1128 PetscCall(KSPDestroy(&user->solver)); 1129 PetscCall(MatDestroy(&user->Js)); 1130 PetscCall(MatDestroy(&user->Jd)); 1131 PetscCall(MatDestroy(&user->JsBlockPrec)); 1132 PetscCall(MatDestroy(&user->JsInv)); 1133 PetscCall(MatDestroy(&user->JsBlock)); 1134 PetscCall(MatDestroy(&user->Divxy[0])); 1135 PetscCall(MatDestroy(&user->Divxy[1])); 1136 PetscCall(MatDestroy(&user->Gradxy[0])); 1137 PetscCall(MatDestroy(&user->Gradxy[1])); 1138 PetscCall(MatDestroy(&user->M)); 1139 for (i = 0; i < user->nt; i++) { 1140 PetscCall(MatDestroy(&user->C[i])); 1141 PetscCall(MatDestroy(&user->Cwork[i])); 1142 } 1143 PetscCall(PetscFree(user->C)); 1144 PetscCall(PetscFree(user->Cwork)); 1145 PetscCall(VecDestroy(&user->u)); 1146 PetscCall(VecDestroy(&user->uwork)); 1147 PetscCall(VecDestroy(&user->vwork)); 1148 PetscCall(VecDestroy(&user->utrue)); 1149 PetscCall(VecDestroy(&user->y)); 1150 PetscCall(VecDestroy(&user->ywork)); 1151 PetscCall(VecDestroy(&user->ytrue)); 1152 PetscCall(VecDestroyVecs(user->nt, &user->yi)); 1153 PetscCall(VecDestroyVecs(user->nt, &user->yiwork)); 1154 PetscCall(VecDestroyVecs(user->nt, &user->ziwork)); 1155 PetscCall(VecDestroyVecs(user->nt, &user->uxi)); 1156 PetscCall(VecDestroyVecs(user->nt, &user->uyi)); 1157 PetscCall(VecDestroyVecs(user->nt, &user->uxiwork)); 1158 PetscCall(VecDestroyVecs(user->nt, &user->uyiwork)); 1159 PetscCall(VecDestroyVecs(user->nt, &user->ui)); 1160 PetscCall(VecDestroyVecs(user->nt, &user->uiwork)); 1161 PetscCall(VecDestroy(&user->c)); 1162 PetscCall(VecDestroy(&user->cwork)); 1163 PetscCall(VecDestroy(&user->ur)); 1164 PetscCall(VecDestroy(&user->q)); 1165 PetscCall(VecDestroy(&user->d)); 1166 PetscCall(VecDestroy(&user->dwork)); 1167 PetscCall(VecDestroy(&user->lwork)); 1168 PetscCall(VecDestroy(&user->js_diag)); 1169 PetscCall(ISDestroy(&user->s_is)); 1170 PetscCall(ISDestroy(&user->d_is)); 1171 PetscCall(VecScatterDestroy(&user->state_scatter)); 1172 PetscCall(VecScatterDestroy(&user->design_scatter)); 1173 for (i = 0; i < user->nt; i++) { 1174 PetscCall(VecScatterDestroy(&user->uxi_scatter[i])); 1175 PetscCall(VecScatterDestroy(&user->uyi_scatter[i])); 1176 PetscCall(VecScatterDestroy(&user->ux_scatter[i])); 1177 PetscCall(VecScatterDestroy(&user->uy_scatter[i])); 1178 PetscCall(VecScatterDestroy(&user->ui_scatter[i])); 1179 PetscCall(VecScatterDestroy(&user->yi_scatter[i])); 1180 } 1181 PetscCall(PetscFree(user->uxi_scatter)); 1182 PetscCall(PetscFree(user->uyi_scatter)); 1183 PetscCall(PetscFree(user->ux_scatter)); 1184 PetscCall(PetscFree(user->uy_scatter)); 1185 PetscCall(PetscFree(user->ui_scatter)); 1186 PetscCall(PetscFree(user->yi_scatter)); 1187 PetscFunctionReturn(PETSC_SUCCESS); 1188 } 1189 1190 PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr) 1191 { 1192 Vec X; 1193 PetscReal unorm, ynorm; 1194 AppCtx *user = (AppCtx *)ptr; 1195 1196 PetscFunctionBegin; 1197 PetscCall(TaoGetSolution(tao, &X)); 1198 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 1199 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue)); 1200 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); 1201 PetscCall(VecNorm(user->uwork, NORM_2, &unorm)); 1202 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm)); 1203 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm)); 1204 PetscFunctionReturn(PETSC_SUCCESS); 1205 } 1206 1207 /*TEST 1208 1209 build: 1210 requires: !complex 1211 1212 test: 1213 requires: !single 1214 args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5 1215 1216 test: 1217 suffix: guess_pod 1218 requires: !single 1219 args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5 1220 1221 TEST*/ 1222