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