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