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