1 #include <petsc/private/taoimpl.h> 2 3 /*T 4 Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares 5 Routines: TaoCreate(); 6 Routines: TaoSetType(); 7 Routines: TaoSetInitialVector(); 8 Routines: TaoSetObjectiveRoutine(); 9 Routines: TaoSetGradientRoutine(); 10 Routines: TaoSetConstraintsRoutine(); 11 Routines: TaoSetJacobianStateRoutine(); 12 Routines: TaoSetJacobianDesignRoutine(); 13 Routines: TaoSetStateDesignIS(); 14 Routines: TaoSetFromOptions(); 15 Routines: TaoSolve(); 16 Routines: TaoDestroy(); 17 Processors: n 18 T*/ 19 20 21 22 typedef struct { 23 PetscInt n; /* Number of total variables */ 24 PetscInt m; /* Number of constraints */ 25 PetscInt nstate; 26 PetscInt ndesign; 27 PetscInt mx; /* grid points in each direction */ 28 PetscInt ns; /* Number of data samples (1<=ns<=8) 29 Currently only ns=1 is supported */ 30 PetscInt ndata; /* Number of data points per sample */ 31 IS s_is; 32 IS d_is; 33 34 VecScatter state_scatter; 35 VecScatter design_scatter; 36 VecScatter *yi_scatter, *di_scatter; 37 Vec suby,subq,subd; 38 Mat Js,Jd,JsPrec,JsInv,JsBlock; 39 40 PetscReal alpha; /* Regularization parameter */ 41 PetscReal beta; /* Weight attributed to ||u||^2 in regularization functional */ 42 PetscReal noise; /* Amount of noise to add to data */ 43 PetscReal *ones; 44 Mat Q; 45 Mat MQ; 46 Mat L; 47 48 Mat Grad; 49 Mat Av,Avwork; 50 Mat Div, Divwork; 51 Mat DSG; 52 Mat Diag,Ones; 53 54 55 Vec q; 56 Vec ur; /* reference */ 57 58 Vec d; 59 Vec dwork; 60 61 Vec x; /* super vec of y,u */ 62 63 Vec y; /* state variables */ 64 Vec ywork; 65 66 Vec ytrue; 67 68 Vec u; /* design variables */ 69 Vec uwork; 70 71 Vec utrue; 72 73 Vec js_diag; 74 75 Vec c; /* constraint vector */ 76 Vec cwork; 77 78 Vec lwork; 79 Vec S; 80 Vec Swork,Twork,Sdiag,Ywork; 81 Vec Av_u; 82 83 KSP solver; 84 PC prec; 85 86 PetscReal tola,tolb,tolc,told; 87 PetscInt ksp_its; 88 PetscInt ksp_its_initial; 89 PetscLogStage stages[10]; 90 PetscBool use_ptap; 91 PetscBool use_lrc; 92 } AppCtx; 93 94 PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*); 95 PetscErrorCode FormGradient(Tao, Vec, Vec, void*); 96 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*); 97 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*); 98 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*); 99 PetscErrorCode FormConstraints(Tao, Vec, Vec, void*); 100 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*); 101 PetscErrorCode Gather(Vec, Vec, VecScatter, Vec, VecScatter); 102 PetscErrorCode Scatter(Vec, Vec, VecScatter, Vec, VecScatter); 103 PetscErrorCode EllipticInitialize(AppCtx*); 104 PetscErrorCode EllipticDestroy(AppCtx*); 105 PetscErrorCode EllipticMonitor(Tao, void*); 106 107 PetscErrorCode StateBlockMatMult(Mat,Vec,Vec); 108 PetscErrorCode StateMatMult(Mat,Vec,Vec); 109 110 PetscErrorCode StateInvMatMult(Mat,Vec,Vec); 111 PetscErrorCode DesignMatMult(Mat,Vec,Vec); 112 PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec); 113 114 PetscErrorCode QMatMult(Mat,Vec,Vec); 115 PetscErrorCode QMatMultTranspose(Mat,Vec,Vec); 116 117 static char help[]=""; 118 119 int main(int argc, char **argv) 120 { 121 PetscErrorCode ierr; 122 Vec x0; 123 Tao tao; 124 AppCtx user; 125 PetscInt ntests = 1; 126 PetscInt i; 127 128 ierr = PetscInitialize(&argc, &argv, (char*)0,help);if (ierr) return ierr; 129 user.mx = 8; 130 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);CHKERRQ(ierr); 131 ierr = PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);CHKERRQ(ierr); 132 user.ns = 6; 133 ierr = PetscOptionsInt("-ns","Number of data samples (1<=ns<=8)","",user.ns,&user.ns,NULL);CHKERRQ(ierr); 134 user.ndata = 64; 135 ierr = PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);CHKERRQ(ierr); 136 user.alpha = 0.1; 137 ierr = PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);CHKERRQ(ierr); 138 user.beta = 0.00001; 139 ierr = PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL);CHKERRQ(ierr); 140 user.noise = 0.01; 141 ierr = PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL);CHKERRQ(ierr); 142 143 user.use_ptap = PETSC_FALSE; 144 ierr = PetscOptionsBool("-use_ptap","Use ptap matrix for DSG","",user.use_ptap,&user.use_ptap,NULL);CHKERRQ(ierr); 145 user.use_lrc = PETSC_FALSE; 146 ierr = PetscOptionsBool("-use_lrc","Use lrc matrix for Js","",user.use_lrc,&user.use_lrc,NULL);CHKERRQ(ierr); 147 user.m = user.ns*user.mx*user.mx*user.mx; /* number of constraints */ 148 user.nstate = user.m; 149 user.ndesign = user.mx*user.mx*user.mx; 150 user.n = user.nstate + user.ndesign; /* number of variables */ 151 152 /* Create TAO solver and set desired solution method */ 153 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 154 ierr = TaoSetType(tao,TAOLCL);CHKERRQ(ierr); 155 156 /* Set up initial vectors and matrices */ 157 ierr = EllipticInitialize(&user);CHKERRQ(ierr); 158 159 ierr = Gather(user.x,user.y,user.state_scatter,user.u,user.design_scatter);CHKERRQ(ierr); 160 ierr = VecDuplicate(user.x,&x0);CHKERRQ(ierr); 161 ierr = VecCopy(user.x,x0);CHKERRQ(ierr); 162 163 /* Set solution vector with an initial guess */ 164 ierr = TaoSetInitialVector(tao,user.x);CHKERRQ(ierr); 165 ierr = TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);CHKERRQ(ierr); 166 ierr = TaoSetGradientRoutine(tao, FormGradient, (void *)&user);CHKERRQ(ierr); 167 ierr = TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);CHKERRQ(ierr); 168 169 ierr = TaoSetJacobianStateRoutine(tao, user.Js, NULL, user.JsInv, FormJacobianState, (void *)&user);CHKERRQ(ierr); 170 ierr = TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);CHKERRQ(ierr); 171 172 ierr = TaoSetStateDesignIS(tao,user.s_is,user.d_is);CHKERRQ(ierr); 173 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 174 175 /* SOLVE THE APPLICATION */ 176 ierr = PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);CHKERRQ(ierr); 177 ierr = PetscOptionsEnd();CHKERRQ(ierr); 178 ierr = PetscLogStageRegister("Trials",&user.stages[1]);CHKERRQ(ierr); 179 ierr = PetscLogStagePush(user.stages[1]);CHKERRQ(ierr); 180 for (i=0; i<ntests; i++){ 181 ierr = TaoSolve(tao);CHKERRQ(ierr); 182 ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its);CHKERRQ(ierr); 183 ierr = VecCopy(x0,user.x);CHKERRQ(ierr); 184 } 185 ierr = PetscLogStagePop();CHKERRQ(ierr); 186 ierr = PetscBarrier((PetscObject)user.x);CHKERRQ(ierr); 187 ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");CHKERRQ(ierr); 188 ierr = PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);CHKERRQ(ierr); 189 190 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 191 ierr = VecDestroy(&x0);CHKERRQ(ierr); 192 ierr = EllipticDestroy(&user);CHKERRQ(ierr); 193 ierr = PetscFinalize(); 194 return ierr; 195 } 196 /* ------------------------------------------------------------------- */ 197 /* 198 dwork = Qy - d 199 lwork = L*(u-ur) 200 f = 1/2 * (dwork.dwork + alpha*lwork.lwork) 201 */ 202 PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr) 203 { 204 PetscErrorCode ierr; 205 PetscReal d1=0,d2=0; 206 AppCtx *user = (AppCtx*)ptr; 207 208 PetscFunctionBegin; 209 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); 210 ierr = MatMult(user->MQ,user->y,user->dwork);CHKERRQ(ierr); 211 ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr); 212 ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr); 213 ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr); 214 ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr); 215 ierr = VecDot(user->lwork,user->lwork,&d2);CHKERRQ(ierr); 216 *f = 0.5 * (d1 + user->alpha*d2); 217 PetscFunctionReturn(0); 218 } 219 220 /* ------------------------------------------------------------------- */ 221 /* 222 state: g_s = Q' *(Qy - d) 223 design: g_d = alpha*L'*L*(u-ur) 224 */ 225 PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr) 226 { 227 PetscErrorCode ierr; 228 AppCtx *user = (AppCtx*)ptr; 229 230 PetscFunctionBegin; 231 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); 232 ierr = MatMult(user->MQ,user->y,user->dwork);CHKERRQ(ierr); 233 ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr); 234 ierr = MatMultTranspose(user->MQ,user->dwork,user->ywork);CHKERRQ(ierr); 235 ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr); 236 ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr); 237 ierr = MatMultTranspose(user->L,user->lwork,user->uwork);CHKERRQ(ierr); 238 ierr = VecScale(user->uwork, user->alpha);CHKERRQ(ierr); 239 ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) 244 { 245 PetscErrorCode ierr; 246 PetscReal d1,d2; 247 AppCtx *user = (AppCtx*)ptr; 248 249 PetscFunctionBegin; 250 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); 251 ierr = MatMult(user->MQ,user->y,user->dwork);CHKERRQ(ierr); 252 ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr); 253 ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr); 254 ierr = MatMultTranspose(user->MQ,user->dwork,user->ywork);CHKERRQ(ierr); 255 256 ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr); 257 ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr); 258 ierr = VecDot(user->lwork,user->lwork,&d2);CHKERRQ(ierr); 259 ierr = MatMultTranspose(user->L,user->lwork,user->uwork);CHKERRQ(ierr); 260 ierr = VecScale(user->uwork, user->alpha);CHKERRQ(ierr); 261 *f = 0.5 * (d1 + user->alpha*d2); 262 ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 /* ------------------------------------------------------------------- */ 267 /* A 268 MatShell object 269 */ 270 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr) 271 { 272 PetscErrorCode ierr; 273 AppCtx *user = (AppCtx*)ptr; 274 275 PetscFunctionBegin; 276 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); 277 /* DSG = Div * (1/Av_u) * Grad */ 278 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); 279 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); 280 ierr = VecExp(user->uwork);CHKERRQ(ierr); 281 ierr = MatMult(user->Av,user->uwork,user->Av_u);CHKERRQ(ierr); 282 ierr = VecCopy(user->Av_u,user->Swork);CHKERRQ(ierr); 283 ierr = VecReciprocal(user->Swork);CHKERRQ(ierr); 284 if (user->use_ptap) { 285 ierr = MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);CHKERRQ(ierr); 286 ierr = MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);CHKERRQ(ierr); 287 } else { 288 ierr = MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 289 ierr = MatDiagonalScale(user->Divwork,NULL,user->Swork);CHKERRQ(ierr); 290 ierr = MatProductNumeric(user->DSG);CHKERRQ(ierr); 291 } 292 PetscFunctionReturn(0); 293 } 294 /* ------------------------------------------------------------------- */ 295 /* B */ 296 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr) 297 { 298 PetscErrorCode ierr; 299 AppCtx *user = (AppCtx*)ptr; 300 301 PetscFunctionBegin; 302 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 PetscErrorCode StateBlockMatMult(Mat J_shell, Vec X, Vec Y) 307 { 308 PetscErrorCode ierr; 309 PetscReal sum; 310 AppCtx *user; 311 312 PetscFunctionBegin; 313 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 314 ierr = MatMult(user->DSG,X,Y);CHKERRQ(ierr); 315 ierr = VecSum(X,&sum);CHKERRQ(ierr); 316 sum /= user->ndesign; 317 ierr = VecShift(Y,sum);CHKERRQ(ierr); 318 PetscFunctionReturn(0); 319 } 320 321 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y) 322 { 323 PetscErrorCode ierr; 324 PetscInt i; 325 AppCtx *user; 326 327 PetscFunctionBegin; 328 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 329 if (user->ns == 1) { 330 ierr = MatMult(user->JsBlock,X,Y);CHKERRQ(ierr); 331 } else { 332 for (i=0;i<user->ns;i++) { 333 ierr = Scatter(X,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 334 ierr = Scatter(Y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 335 ierr = MatMult(user->JsBlock,user->subq,user->suby);CHKERRQ(ierr); 336 ierr = Gather(Y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 337 } 338 } 339 PetscFunctionReturn(0); 340 } 341 342 PetscErrorCode StateInvMatMult(Mat J_shell, Vec X, Vec Y) 343 { 344 PetscErrorCode ierr; 345 PetscInt its,i; 346 AppCtx *user; 347 348 PetscFunctionBegin; 349 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 350 ierr = KSPSetOperators(user->solver,user->JsBlock,user->DSG);CHKERRQ(ierr); 351 if (Y == user->ytrue) { 352 /* First solve is done using true solution to set up problem */ 353 ierr = KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 354 } else { 355 ierr = KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 356 } 357 if (user->ns == 1) { 358 ierr = KSPSolve(user->solver,X,Y);CHKERRQ(ierr); 359 ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr); 360 user->ksp_its+=its; 361 } else { 362 for (i=0;i<user->ns;i++) { 363 ierr = Scatter(X,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 364 ierr = Scatter(Y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 365 ierr = KSPSolve(user->solver,user->subq,user->suby);CHKERRQ(ierr); 366 ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr); 367 user->ksp_its+=its; 368 ierr = Gather(Y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 369 } 370 } 371 PetscFunctionReturn(0); 372 } 373 PetscErrorCode QMatMult(Mat J_shell, Vec X, Vec Y) 374 { 375 PetscErrorCode ierr; 376 AppCtx *user; 377 PetscInt i; 378 379 PetscFunctionBegin; 380 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 381 if (user->ns == 1) { 382 ierr = MatMult(user->Q,X,Y);CHKERRQ(ierr); 383 } else { 384 for (i=0;i<user->ns;i++) { 385 ierr = Scatter(X,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 386 ierr = Scatter(Y,user->subd,user->di_scatter[i],0,0);CHKERRQ(ierr); 387 ierr = MatMult(user->Q,user->subq,user->subd);CHKERRQ(ierr); 388 ierr = Gather(Y,user->subd,user->di_scatter[i],0,0);CHKERRQ(ierr); 389 } 390 } 391 PetscFunctionReturn(0); 392 } 393 394 PetscErrorCode QMatMultTranspose(Mat J_shell, Vec X, Vec Y) 395 { 396 PetscErrorCode ierr; 397 AppCtx *user; 398 PetscInt i; 399 400 PetscFunctionBegin; 401 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 402 if (user->ns == 1) { 403 ierr = MatMultTranspose(user->Q,X,Y);CHKERRQ(ierr); 404 } else { 405 for (i=0;i<user->ns;i++) { 406 ierr = Scatter(X,user->subd,user->di_scatter[i],0,0);CHKERRQ(ierr); 407 ierr = Scatter(Y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 408 ierr = MatMultTranspose(user->Q,user->subd,user->suby);CHKERRQ(ierr); 409 ierr = Gather(Y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 410 } 411 } 412 PetscFunctionReturn(0); 413 } 414 415 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y) 416 { 417 PetscErrorCode ierr; 418 PetscInt i; 419 AppCtx *user; 420 421 PetscFunctionBegin; 422 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 423 424 /* sdiag(1./v) */ 425 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); 426 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); 427 ierr = VecExp(user->uwork);CHKERRQ(ierr); 428 429 /* sdiag(1./((Av*(1./v)).^2)) */ 430 ierr = MatMult(user->Av,user->uwork,user->Swork);CHKERRQ(ierr); 431 ierr = VecPointwiseMult(user->Swork,user->Swork,user->Swork);CHKERRQ(ierr); 432 ierr = VecReciprocal(user->Swork);CHKERRQ(ierr); 433 434 /* (Av * (sdiag(1./v) * b)) */ 435 ierr = VecPointwiseMult(user->uwork,user->uwork,X);CHKERRQ(ierr); 436 ierr = MatMult(user->Av,user->uwork,user->Twork);CHKERRQ(ierr); 437 438 /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */ 439 ierr = VecPointwiseMult(user->Swork,user->Twork,user->Swork);CHKERRQ(ierr); 440 441 if (user->ns == 1) { 442 /* (sdiag(Grad*y(:,i)) */ 443 ierr = MatMult(user->Grad,user->y,user->Twork);CHKERRQ(ierr); 444 445 /* Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */ 446 ierr = VecPointwiseMult(user->Swork,user->Twork,user->Swork);CHKERRQ(ierr); 447 ierr = MatMultTranspose(user->Grad,user->Swork,Y);CHKERRQ(ierr); 448 } else { 449 for (i=0;i<user->ns;i++) { 450 ierr = Scatter(user->y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 451 ierr = Scatter(Y,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 452 453 ierr = MatMult(user->Grad,user->suby,user->Twork);CHKERRQ(ierr); 454 ierr = VecPointwiseMult(user->Twork,user->Twork,user->Swork);CHKERRQ(ierr); 455 ierr = MatMultTranspose(user->Grad,user->Twork,user->subq);CHKERRQ(ierr); 456 ierr = Gather(user->y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 457 ierr = Gather(Y,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 458 } 459 } 460 PetscFunctionReturn(0); 461 } 462 463 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 464 { 465 PetscErrorCode ierr; 466 PetscInt i; 467 AppCtx *user; 468 469 PetscFunctionBegin; 470 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); 471 ierr = VecZeroEntries(Y);CHKERRQ(ierr); 472 473 /* Sdiag = 1./((Av*(1./v)).^2) */ 474 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); 475 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); 476 ierr = VecExp(user->uwork);CHKERRQ(ierr); 477 ierr = MatMult(user->Av,user->uwork,user->Swork);CHKERRQ(ierr); 478 ierr = VecPointwiseMult(user->Sdiag,user->Swork,user->Swork);CHKERRQ(ierr); 479 ierr = VecReciprocal(user->Sdiag);CHKERRQ(ierr); 480 481 for (i=0;i<user->ns;i++) { 482 ierr = Scatter(X,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 483 ierr = Scatter(user->y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 484 485 /* Swork = (Div' * b(:,i)) */ 486 ierr = MatMult(user->Grad,user->subq,user->Swork);CHKERRQ(ierr); 487 488 /* Twork = Grad*y(:,i) */ 489 ierr = MatMult(user->Grad,user->suby,user->Twork);CHKERRQ(ierr); 490 491 /* Twork = sdiag(Twork) * Swork */ 492 ierr = VecPointwiseMult(user->Twork,user->Swork,user->Twork);CHKERRQ(ierr); 493 494 495 /* Swork = pointwisemult(Sdiag,Twork) */ 496 ierr = VecPointwiseMult(user->Swork,user->Twork,user->Sdiag);CHKERRQ(ierr); 497 498 /* Ywork = Av' * Swork */ 499 ierr = MatMultTranspose(user->Av,user->Swork,user->Ywork);CHKERRQ(ierr); 500 501 /* Ywork = pointwisemult(uwork,Ywork) */ 502 ierr = VecPointwiseMult(user->Ywork,user->uwork,user->Ywork);CHKERRQ(ierr); 503 ierr = VecAXPY(Y,1.0,user->Ywork);CHKERRQ(ierr); 504 ierr = Gather(user->y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 505 } 506 PetscFunctionReturn(0); 507 } 508 509 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr) 510 { 511 /* C=Ay - q A = Div * Sigma * Grad + hx*hx*hx*ones(n,n) */ 512 PetscErrorCode ierr; 513 PetscReal sum; 514 PetscInt i; 515 AppCtx *user = (AppCtx*)ptr; 516 517 PetscFunctionBegin; 518 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); 519 if (user->ns == 1) { 520 ierr = MatMult(user->Grad,user->y,user->Swork);CHKERRQ(ierr); 521 ierr = VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);CHKERRQ(ierr); 522 ierr = MatMultTranspose(user->Grad,user->Swork,C);CHKERRQ(ierr); 523 ierr = VecSum(user->y,&sum);CHKERRQ(ierr); 524 sum /= user->ndesign; 525 ierr = VecShift(C,sum);CHKERRQ(ierr); 526 } else { 527 for (i=0;i<user->ns;i++) { 528 ierr = Scatter(user->y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 529 ierr = Scatter(C,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 530 ierr = MatMult(user->Grad,user->suby,user->Swork);CHKERRQ(ierr); 531 ierr = VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);CHKERRQ(ierr); 532 ierr = MatMultTranspose(user->Grad,user->Swork,user->subq);CHKERRQ(ierr); 533 534 ierr = VecSum(user->suby,&sum);CHKERRQ(ierr); 535 sum /= user->ndesign; 536 ierr = VecShift(user->subq,sum);CHKERRQ(ierr); 537 538 ierr = Gather(user->y,user->suby,user->yi_scatter[i],0,0);CHKERRQ(ierr); 539 ierr = Gather(C,user->subq,user->yi_scatter[i],0,0);CHKERRQ(ierr); 540 } 541 } 542 ierr = VecAXPY(C,-1.0,user->q);CHKERRQ(ierr); 543 PetscFunctionReturn(0); 544 } 545 546 PetscErrorCode Scatter(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2) 547 { 548 PetscErrorCode ierr; 549 550 PetscFunctionBegin; 551 ierr = VecScatterBegin(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 552 ierr = VecScatterEnd(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 553 if (sub2) { 554 ierr = VecScatterBegin(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 555 ierr = VecScatterEnd(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 556 } 557 PetscFunctionReturn(0); 558 } 559 560 PetscErrorCode Gather(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2) 561 { 562 PetscErrorCode ierr; 563 564 PetscFunctionBegin; 565 ierr = VecScatterBegin(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 566 ierr = VecScatterEnd(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 567 if (sub2) { 568 ierr = VecScatterBegin(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 569 ierr = VecScatterEnd(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 570 } 571 PetscFunctionReturn(0); 572 } 573 574 PetscErrorCode EllipticInitialize(AppCtx *user) 575 { 576 PetscErrorCode ierr; 577 PetscInt m,n,i,j,k,l,linear_index,is,js,ks,ls,istart,iend,iblock; 578 Vec XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork; 579 PetscReal *x,*y,*z; 580 PetscReal h,meanut; 581 PetscScalar hinv,neg_hinv,half = 0.5,sqrt_beta; 582 PetscInt im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz; 583 IS is_alldesign,is_allstate; 584 IS is_from_d; 585 IS is_from_y; 586 PetscInt lo,hi,hi2,lo2,ysubnlocal,dsubnlocal; 587 const PetscInt *ranges, *subranges; 588 PetscMPIInt size; 589 PetscReal xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz; 590 PetscScalar v,vx,vy,vz; 591 PetscInt offset,subindex,subvec,nrank,kk; 592 593 PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160, 594 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774, 595 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997, 0.5198, 0.8326, 596 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728, 597 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273, 598 0.5724, 0.4331, 0.5136, 0.3547, 0.4413, 0.2602, 0.5698, 0.7278, 599 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594, 600 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756}; 601 602 PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256, 603 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792, 604 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437, 0.5938, 0.6137, 605 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985, 606 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288, 607 0.5761, 0.1129, 0.3841, 0.6150, 0.6904, 0.6686, 0.1361, 0.4601, 608 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480, 609 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658}; 610 611 PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295, 612 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819, 613 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236, 0.8800, 0.2939, 614 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712, 615 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078, 616 0.1987, 0.2297, 0.4321, 0.8115, 0.4915, 0.7764, 0.4657, 0.4627, 617 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313, 618 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711}; 619 620 PetscFunctionBegin; 621 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 622 ierr = PetscLogStageRegister("Elliptic Setup",&user->stages[0]);CHKERRQ(ierr); 623 ierr = PetscLogStagePush(user->stages[0]);CHKERRQ(ierr); 624 625 /* Create u,y,c,x */ 626 ierr = VecCreate(PETSC_COMM_WORLD,&user->u);CHKERRQ(ierr); 627 ierr = VecCreate(PETSC_COMM_WORLD,&user->y);CHKERRQ(ierr); 628 ierr = VecCreate(PETSC_COMM_WORLD,&user->c);CHKERRQ(ierr); 629 ierr = VecSetSizes(user->u,PETSC_DECIDE,user->ndesign);CHKERRQ(ierr); 630 ierr = VecSetFromOptions(user->u);CHKERRQ(ierr); 631 ierr = VecGetLocalSize(user->u,&ysubnlocal);CHKERRQ(ierr); 632 ierr = VecSetSizes(user->y,ysubnlocal*user->ns,user->nstate);CHKERRQ(ierr); 633 ierr = VecSetSizes(user->c,ysubnlocal*user->ns,user->m);CHKERRQ(ierr); 634 ierr = VecSetFromOptions(user->y);CHKERRQ(ierr); 635 ierr = VecSetFromOptions(user->c);CHKERRQ(ierr); 636 637 /* 638 ******************************* 639 Create scatters for x <-> y,u 640 ******************************* 641 642 If the state vector y and design vector u are partitioned as 643 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors), 644 then the solution vector x is organized as 645 [y_1; u_1; y_2; u_2; ...; y_np; u_np]. 646 The index sets user->s_is and user->d_is correspond to the indices of the 647 state and design variables owned by the current processor. 648 */ 649 ierr = VecCreate(PETSC_COMM_WORLD,&user->x);CHKERRQ(ierr); 650 651 ierr = VecGetOwnershipRange(user->y,&lo,&hi);CHKERRQ(ierr); 652 ierr = VecGetOwnershipRange(user->u,&lo2,&hi2);CHKERRQ(ierr); 653 654 ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);CHKERRQ(ierr); 655 ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user->s_is);CHKERRQ(ierr); 656 ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);CHKERRQ(ierr); 657 ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user->d_is);CHKERRQ(ierr); 658 659 ierr = VecSetSizes(user->x,hi-lo+hi2-lo2,user->n);CHKERRQ(ierr); 660 ierr = VecSetFromOptions(user->x);CHKERRQ(ierr); 661 662 ierr = VecScatterCreate(user->x,user->s_is,user->y,is_allstate,&user->state_scatter);CHKERRQ(ierr); 663 ierr = VecScatterCreate(user->x,user->d_is,user->u,is_alldesign,&user->design_scatter);CHKERRQ(ierr); 664 ierr = ISDestroy(&is_alldesign);CHKERRQ(ierr); 665 ierr = ISDestroy(&is_allstate);CHKERRQ(ierr); 666 /* 667 ******************************* 668 Create scatter from y to y_1,y_2,...,y_ns 669 ******************************* 670 */ 671 ierr = PetscMalloc1(user->ns,&user->yi_scatter);CHKERRQ(ierr); 672 ierr = VecDuplicate(user->u,&user->suby);CHKERRQ(ierr); 673 ierr = VecDuplicate(user->u,&user->subq);CHKERRQ(ierr); 674 675 ierr = VecGetOwnershipRange(user->y,&lo2,&hi2);CHKERRQ(ierr); 676 istart = 0; 677 for (i=0; i<user->ns; i++){ 678 ierr = VecGetOwnershipRange(user->suby,&lo,&hi);CHKERRQ(ierr); 679 ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y);CHKERRQ(ierr); 680 ierr = VecScatterCreate(user->y,is_from_y,user->suby,NULL,&user->yi_scatter[i]);CHKERRQ(ierr); 681 istart = istart + hi-lo; 682 ierr = ISDestroy(&is_from_y);CHKERRQ(ierr); 683 } 684 /* 685 ******************************* 686 Create scatter from d to d_1,d_2,...,d_ns 687 ******************************* 688 */ 689 ierr = VecCreate(PETSC_COMM_WORLD,&user->subd);CHKERRQ(ierr); 690 ierr = VecSetSizes(user->subd,PETSC_DECIDE,user->ndata);CHKERRQ(ierr); 691 ierr = VecSetFromOptions(user->subd);CHKERRQ(ierr); 692 ierr = VecCreate(PETSC_COMM_WORLD,&user->d);CHKERRQ(ierr); 693 ierr = VecGetLocalSize(user->subd,&dsubnlocal);CHKERRQ(ierr); 694 ierr = VecSetSizes(user->d,dsubnlocal*user->ns,user->ndata*user->ns);CHKERRQ(ierr); 695 ierr = VecSetFromOptions(user->d);CHKERRQ(ierr); 696 ierr = PetscMalloc1(user->ns,&user->di_scatter);CHKERRQ(ierr); 697 698 ierr = VecGetOwnershipRange(user->d,&lo2,&hi2);CHKERRQ(ierr); 699 istart = 0; 700 for (i=0; i<user->ns; i++){ 701 ierr = VecGetOwnershipRange(user->subd,&lo,&hi);CHKERRQ(ierr); 702 ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d);CHKERRQ(ierr); 703 ierr = VecScatterCreate(user->d,is_from_d,user->subd,NULL,&user->di_scatter[i]);CHKERRQ(ierr); 704 istart = istart + hi-lo; 705 ierr = ISDestroy(&is_from_d);CHKERRQ(ierr); 706 } 707 708 ierr = PetscMalloc1(user->mx,&x);CHKERRQ(ierr); 709 ierr = PetscMalloc1(user->mx,&y);CHKERRQ(ierr); 710 ierr = PetscMalloc1(user->mx,&z);CHKERRQ(ierr); 711 712 user->ksp_its = 0; 713 user->ksp_its_initial = 0; 714 715 n = user->mx * user->mx * user->mx; 716 m = 3 * user->mx * user->mx * (user->mx-1); 717 sqrt_beta = PetscSqrtScalar(user->beta); 718 719 ierr = VecCreate(PETSC_COMM_WORLD,&XX);CHKERRQ(ierr); 720 ierr = VecCreate(PETSC_COMM_WORLD,&user->q);CHKERRQ(ierr); 721 ierr = VecSetSizes(XX,ysubnlocal,n);CHKERRQ(ierr); 722 ierr = VecSetSizes(user->q,ysubnlocal*user->ns,user->m);CHKERRQ(ierr); 723 ierr = VecSetFromOptions(XX);CHKERRQ(ierr); 724 ierr = VecSetFromOptions(user->q);CHKERRQ(ierr); 725 726 ierr = VecDuplicate(XX,&YY);CHKERRQ(ierr); 727 ierr = VecDuplicate(XX,&ZZ);CHKERRQ(ierr); 728 ierr = VecDuplicate(XX,&XXwork);CHKERRQ(ierr); 729 ierr = VecDuplicate(XX,&YYwork);CHKERRQ(ierr); 730 ierr = VecDuplicate(XX,&ZZwork);CHKERRQ(ierr); 731 ierr = VecDuplicate(XX,&UTwork);CHKERRQ(ierr); 732 ierr = VecDuplicate(XX,&user->utrue);CHKERRQ(ierr); 733 734 /* map for striding q */ 735 ierr = VecGetOwnershipRanges(user->q,&ranges);CHKERRQ(ierr); 736 ierr = VecGetOwnershipRanges(user->u,&subranges);CHKERRQ(ierr); 737 738 ierr = VecGetOwnershipRange(user->q,&lo2,&hi2);CHKERRQ(ierr); 739 ierr = VecGetOwnershipRange(user->u,&lo,&hi);CHKERRQ(ierr); 740 /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */ 741 h = 1.0/user->mx; 742 hinv = user->mx; 743 neg_hinv = -hinv; 744 745 ierr = VecGetOwnershipRange(XX,&istart,&iend);CHKERRQ(ierr); 746 for (linear_index=istart; linear_index<iend; linear_index++){ 747 i = linear_index % user->mx; 748 j = ((linear_index-i)/user->mx) % user->mx; 749 k = ((linear_index-i)/user->mx-j) / user->mx; 750 vx = h*(i+0.5); 751 vy = h*(j+0.5); 752 vz = h*(k+0.5); 753 ierr = VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);CHKERRQ(ierr); 754 ierr = VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);CHKERRQ(ierr); 755 ierr = VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES);CHKERRQ(ierr); 756 for (is=0; is<2; is++){ 757 for (js=0; js<2; js++){ 758 for (ks=0; ks<2; ks++){ 759 ls = is*4 + js*2 + ks; 760 if (ls<user->ns){ 761 l =ls*n + linear_index; 762 /* remap */ 763 subindex = l%n; 764 subvec = l/n; 765 nrank=0; 766 while (subindex >= subranges[nrank+1]) nrank++; 767 offset = subindex - subranges[nrank]; 768 istart=0; 769 for (kk=0;kk<nrank;kk++) istart+=user->ns*(subranges[kk+1]-subranges[kk]); 770 istart += (subranges[nrank+1]-subranges[nrank])*subvec; 771 l = istart+offset; 772 v = 100*PetscSinScalar(2*PETSC_PI*(vx+0.25*is))*PetscSinScalar(2*PETSC_PI*(vy+0.25*js))*PetscSinScalar(2*PETSC_PI*(vz+0.25*ks)); 773 ierr = VecSetValues(user->q,1,&l,&v,INSERT_VALUES);CHKERRQ(ierr); 774 } 775 } 776 } 777 } 778 } 779 780 ierr = VecAssemblyBegin(XX);CHKERRQ(ierr); 781 ierr = VecAssemblyEnd(XX);CHKERRQ(ierr); 782 ierr = VecAssemblyBegin(YY);CHKERRQ(ierr); 783 ierr = VecAssemblyEnd(YY);CHKERRQ(ierr); 784 ierr = VecAssemblyBegin(ZZ);CHKERRQ(ierr); 785 ierr = VecAssemblyEnd(ZZ);CHKERRQ(ierr); 786 ierr = VecAssemblyBegin(user->q);CHKERRQ(ierr); 787 ierr = VecAssemblyEnd(user->q);CHKERRQ(ierr); 788 789 /* Compute true parameter function 790 ut = exp(-((x-0.25)^2+(y-0.25)^2+(z-0.25)^2)/0.05) - exp((x-0.75)^2-(y-0.75)^2-(z-0.75))^2/0.05) */ 791 ierr = VecCopy(XX,XXwork);CHKERRQ(ierr); 792 ierr = VecCopy(YY,YYwork);CHKERRQ(ierr); 793 ierr = VecCopy(ZZ,ZZwork);CHKERRQ(ierr); 794 795 ierr = VecShift(XXwork,-0.25);CHKERRQ(ierr); 796 ierr = VecShift(YYwork,-0.25);CHKERRQ(ierr); 797 ierr = VecShift(ZZwork,-0.25);CHKERRQ(ierr); 798 799 ierr = VecPointwiseMult(XXwork,XXwork,XXwork);CHKERRQ(ierr); 800 ierr = VecPointwiseMult(YYwork,YYwork,YYwork);CHKERRQ(ierr); 801 ierr = VecPointwiseMult(ZZwork,ZZwork,ZZwork);CHKERRQ(ierr); 802 803 ierr = VecCopy(XXwork,UTwork);CHKERRQ(ierr); 804 ierr = VecAXPY(UTwork,1.0,YYwork);CHKERRQ(ierr); 805 ierr = VecAXPY(UTwork,1.0,ZZwork);CHKERRQ(ierr); 806 ierr = VecScale(UTwork,-20.0);CHKERRQ(ierr); 807 ierr = VecExp(UTwork);CHKERRQ(ierr); 808 ierr = VecCopy(UTwork,user->utrue);CHKERRQ(ierr); 809 810 ierr = VecCopy(XX,XXwork);CHKERRQ(ierr); 811 ierr = VecCopy(YY,YYwork);CHKERRQ(ierr); 812 ierr = VecCopy(ZZ,ZZwork);CHKERRQ(ierr); 813 814 ierr = VecShift(XXwork,-0.75);CHKERRQ(ierr); 815 ierr = VecShift(YYwork,-0.75);CHKERRQ(ierr); 816 ierr = VecShift(ZZwork,-0.75);CHKERRQ(ierr); 817 818 ierr = VecPointwiseMult(XXwork,XXwork,XXwork);CHKERRQ(ierr); 819 ierr = VecPointwiseMult(YYwork,YYwork,YYwork);CHKERRQ(ierr); 820 ierr = VecPointwiseMult(ZZwork,ZZwork,ZZwork);CHKERRQ(ierr); 821 822 ierr = VecCopy(XXwork,UTwork);CHKERRQ(ierr); 823 ierr = VecAXPY(UTwork,1.0,YYwork);CHKERRQ(ierr); 824 ierr = VecAXPY(UTwork,1.0,ZZwork);CHKERRQ(ierr); 825 ierr = VecScale(UTwork,-20.0);CHKERRQ(ierr); 826 ierr = VecExp(UTwork);CHKERRQ(ierr); 827 828 ierr = VecAXPY(user->utrue,-1.0,UTwork);CHKERRQ(ierr); 829 830 ierr = VecDestroy(&XX);CHKERRQ(ierr); 831 ierr = VecDestroy(&YY);CHKERRQ(ierr); 832 ierr = VecDestroy(&ZZ);CHKERRQ(ierr); 833 ierr = VecDestroy(&XXwork);CHKERRQ(ierr); 834 ierr = VecDestroy(&YYwork);CHKERRQ(ierr); 835 ierr = VecDestroy(&ZZwork);CHKERRQ(ierr); 836 ierr = VecDestroy(&UTwork);CHKERRQ(ierr); 837 838 /* Initial guess and reference model */ 839 ierr = VecDuplicate(user->utrue,&user->ur);CHKERRQ(ierr); 840 ierr = VecSum(user->utrue,&meanut);CHKERRQ(ierr); 841 meanut = meanut / n; 842 ierr = VecSet(user->ur,meanut);CHKERRQ(ierr); 843 ierr = VecCopy(user->ur,user->u);CHKERRQ(ierr); 844 845 /* Generate Grad matrix */ 846 ierr = MatCreate(PETSC_COMM_WORLD,&user->Grad);CHKERRQ(ierr); 847 ierr = MatSetSizes(user->Grad,PETSC_DECIDE,ysubnlocal,m,n);CHKERRQ(ierr); 848 ierr = MatSetFromOptions(user->Grad);CHKERRQ(ierr); 849 ierr = MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL);CHKERRQ(ierr); 850 ierr = MatSeqAIJSetPreallocation(user->Grad,2,NULL);CHKERRQ(ierr); 851 ierr = MatGetOwnershipRange(user->Grad,&istart,&iend);CHKERRQ(ierr); 852 853 for (i=istart; i<iend; i++){ 854 if (i<m/3){ 855 iblock = i / (user->mx-1); 856 j = iblock*user->mx + (i % (user->mx-1)); 857 ierr = MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 858 j = j+1; 859 ierr = MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); 860 } 861 if (i>=m/3 && i<2*m/3){ 862 iblock = (i-m/3) / (user->mx*(user->mx-1)); 863 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 864 ierr = MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 865 j = j + user->mx; 866 ierr = MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); 867 } 868 if (i>=2*m/3){ 869 j = i-2*m/3; 870 ierr = MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 871 j = j + user->mx*user->mx; 872 ierr = MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); 873 } 874 } 875 876 ierr = MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 877 ierr = MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 878 879 /* Generate arithmetic averaging matrix Av */ 880 ierr = MatCreate(PETSC_COMM_WORLD,&user->Av);CHKERRQ(ierr); 881 ierr = MatSetSizes(user->Av,PETSC_DECIDE,ysubnlocal,m,n);CHKERRQ(ierr); 882 ierr = MatSetFromOptions(user->Av);CHKERRQ(ierr); 883 ierr = MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL);CHKERRQ(ierr); 884 ierr = MatSeqAIJSetPreallocation(user->Av,2,NULL);CHKERRQ(ierr); 885 ierr = MatGetOwnershipRange(user->Av,&istart,&iend);CHKERRQ(ierr); 886 887 for (i=istart; i<iend; i++){ 888 if (i<m/3){ 889 iblock = i / (user->mx-1); 890 j = iblock*user->mx + (i % (user->mx-1)); 891 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); 892 j = j+1; 893 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); 894 } 895 if (i>=m/3 && i<2*m/3){ 896 iblock = (i-m/3) / (user->mx*(user->mx-1)); 897 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 898 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); 899 j = j + user->mx; 900 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); 901 } 902 if (i>=2*m/3){ 903 j = i-2*m/3; 904 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); 905 j = j + user->mx*user->mx; 906 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); 907 } 908 } 909 910 ierr = MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 911 ierr = MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 912 913 ierr = MatCreate(PETSC_COMM_WORLD,&user->L);CHKERRQ(ierr); 914 ierr = MatSetSizes(user->L,PETSC_DECIDE,ysubnlocal,m+n,n);CHKERRQ(ierr); 915 ierr = MatSetFromOptions(user->L);CHKERRQ(ierr); 916 ierr = MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL);CHKERRQ(ierr); 917 ierr = MatSeqAIJSetPreallocation(user->L,2,NULL);CHKERRQ(ierr); 918 ierr = MatGetOwnershipRange(user->L,&istart,&iend);CHKERRQ(ierr); 919 920 for (i=istart; i<iend; i++){ 921 if (i<m/3){ 922 iblock = i / (user->mx-1); 923 j = iblock*user->mx + (i % (user->mx-1)); 924 ierr = MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 925 j = j+1; 926 ierr = MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); 927 } 928 if (i>=m/3 && i<2*m/3){ 929 iblock = (i-m/3) / (user->mx*(user->mx-1)); 930 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 931 ierr = MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 932 j = j + user->mx; 933 ierr = MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); 934 } 935 if (i>=2*m/3 && i<m){ 936 j = i-2*m/3; 937 ierr = MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 938 j = j + user->mx*user->mx; 939 ierr = MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); 940 } 941 if (i>=m){ 942 j = i - m; 943 ierr = MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES);CHKERRQ(ierr); 944 } 945 } 946 ierr = MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 947 ierr = MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 948 ierr = MatScale(user->L,PetscPowScalar(h,1.5));CHKERRQ(ierr); 949 950 /* Generate Div matrix */ 951 if (!user->use_ptap) { 952 /* Generate Div matrix */ 953 ierr = MatCreate(PETSC_COMM_WORLD,&user->Div);CHKERRQ(ierr); 954 ierr = MatSetSizes(user->Div,ysubnlocal,PETSC_DECIDE,n,m);CHKERRQ(ierr); 955 ierr = MatSetFromOptions(user->Div);CHKERRQ(ierr); 956 ierr = MatMPIAIJSetPreallocation(user->Div,4,NULL,4,NULL);CHKERRQ(ierr); 957 ierr = MatSeqAIJSetPreallocation(user->Div,6,NULL);CHKERRQ(ierr); 958 ierr = MatGetOwnershipRange(user->Grad,&istart,&iend);CHKERRQ(ierr); 959 960 for (i=istart; i<iend; i++){ 961 if (i<m/3){ 962 iblock = i / (user->mx-1); 963 j = iblock*user->mx + (i % (user->mx-1)); 964 ierr = MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 965 j = j+1; 966 ierr = MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);CHKERRQ(ierr); 967 } 968 if (i>=m/3 && i<2*m/3){ 969 iblock = (i-m/3) / (user->mx*(user->mx-1)); 970 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 971 ierr = MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 972 j = j + user->mx; 973 ierr = MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);CHKERRQ(ierr); 974 } 975 if (i>=2*m/3){ 976 j = i-2*m/3; 977 ierr = MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); 978 j = j + user->mx*user->mx; 979 ierr = MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);CHKERRQ(ierr); 980 } 981 } 982 983 ierr = MatAssemblyBegin(user->Div,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 984 ierr = MatAssemblyEnd(user->Div,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 985 ierr = MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);CHKERRQ(ierr); 986 } else { 987 ierr = MatCreate(PETSC_COMM_WORLD,&user->Diag);CHKERRQ(ierr); 988 ierr = MatSetSizes(user->Diag,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(ierr); 989 ierr = MatSetFromOptions(user->Diag);CHKERRQ(ierr); 990 ierr = MatMPIAIJSetPreallocation(user->Diag,1,NULL,0,NULL);CHKERRQ(ierr); 991 ierr = MatSeqAIJSetPreallocation(user->Diag,1,NULL);CHKERRQ(ierr); 992 } 993 994 /* Build work vectors and matrices */ 995 ierr = VecCreate(PETSC_COMM_WORLD,&user->S);CHKERRQ(ierr); 996 ierr = VecSetSizes(user->S, PETSC_DECIDE, m);CHKERRQ(ierr); 997 ierr = VecSetFromOptions(user->S);CHKERRQ(ierr); 998 999 ierr = VecCreate(PETSC_COMM_WORLD,&user->lwork);CHKERRQ(ierr); 1000 ierr = VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);CHKERRQ(ierr); 1001 ierr = VecSetFromOptions(user->lwork);CHKERRQ(ierr); 1002 1003 ierr = MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork);CHKERRQ(ierr); 1004 1005 ierr = VecDuplicate(user->S,&user->Swork);CHKERRQ(ierr); 1006 ierr = VecDuplicate(user->S,&user->Sdiag);CHKERRQ(ierr); 1007 ierr = VecDuplicate(user->S,&user->Av_u);CHKERRQ(ierr); 1008 ierr = VecDuplicate(user->S,&user->Twork);CHKERRQ(ierr); 1009 ierr = VecDuplicate(user->y,&user->ywork);CHKERRQ(ierr); 1010 ierr = VecDuplicate(user->u,&user->Ywork);CHKERRQ(ierr); 1011 ierr = VecDuplicate(user->u,&user->uwork);CHKERRQ(ierr); 1012 ierr = VecDuplicate(user->u,&user->js_diag);CHKERRQ(ierr); 1013 ierr = VecDuplicate(user->c,&user->cwork);CHKERRQ(ierr); 1014 ierr = VecDuplicate(user->d,&user->dwork);CHKERRQ(ierr); 1015 1016 /* Create a matrix-free shell user->Jd for computing B*x */ 1017 ierr = MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal,user->nstate,user->ndesign,user,&user->Jd);CHKERRQ(ierr); 1018 ierr = MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);CHKERRQ(ierr); 1019 ierr = MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);CHKERRQ(ierr); 1020 1021 /* Compute true state function ytrue given utrue */ 1022 ierr = VecDuplicate(user->y,&user->ytrue);CHKERRQ(ierr); 1023 1024 /* First compute Av_u = Av*exp(-u) */ 1025 ierr = VecSet(user->uwork, 0);CHKERRQ(ierr); 1026 ierr = VecAXPY(user->uwork,-1.0,user->utrue);CHKERRQ(ierr); /* Note: user->utrue */ 1027 ierr = VecExp(user->uwork);CHKERRQ(ierr); 1028 ierr = MatMult(user->Av,user->uwork,user->Av_u);CHKERRQ(ierr); 1029 1030 /* Next form DSG = Div*S*Grad */ 1031 ierr = VecCopy(user->Av_u,user->Swork);CHKERRQ(ierr); 1032 ierr = VecReciprocal(user->Swork);CHKERRQ(ierr); 1033 if (user->use_ptap) { 1034 ierr = MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);CHKERRQ(ierr); 1035 ierr = MatPtAP(user->Diag,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);CHKERRQ(ierr); 1036 } else { 1037 ierr = MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1038 ierr = MatDiagonalScale(user->Divwork,NULL,user->Swork);CHKERRQ(ierr); 1039 1040 ierr = MatMatMult(user->Divwork,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);CHKERRQ(ierr); 1041 } 1042 1043 ierr = MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 1044 ierr = MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); 1045 1046 if (user->use_lrc == PETSC_TRUE) { 1047 v=PetscSqrtReal(1.0 /user->ndesign); 1048 ierr = PetscMalloc1(user->ndesign,&user->ones);CHKERRQ(ierr); 1049 1050 for (i=0;i<user->ndesign;i++) { 1051 user->ones[i]=v; 1052 } 1053 ierr = MatCreateDense(PETSC_COMM_WORLD,ysubnlocal,PETSC_DECIDE,user->ndesign,1,user->ones,&user->Ones);CHKERRQ(ierr); 1054 ierr = MatAssemblyBegin(user->Ones, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1055 ierr = MatAssemblyEnd(user->Ones, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1056 ierr = MatCreateLRC(user->DSG,user->Ones,NULL,user->Ones,&user->JsBlock);CHKERRQ(ierr); 1057 ierr = MatSetUp(user->JsBlock);CHKERRQ(ierr); 1058 } else { 1059 /* Create matrix-free shell user->Js for computing (A + h^3*e*e^T)*x */ 1060 ierr = MatCreateShell(PETSC_COMM_WORLD,ysubnlocal,ysubnlocal,user->ndesign,user->ndesign,user,&user->JsBlock);CHKERRQ(ierr); 1061 ierr = MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateBlockMatMult);CHKERRQ(ierr); 1062 ierr = MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateBlockMatMult);CHKERRQ(ierr); 1063 } 1064 ierr = MatSetOption(user->JsBlock,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 1065 ierr = MatSetOption(user->JsBlock,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); 1066 ierr = MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->Js);CHKERRQ(ierr); 1067 ierr = MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);CHKERRQ(ierr); 1068 ierr = MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMult);CHKERRQ(ierr); 1069 ierr = MatSetOption(user->Js,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 1070 ierr = MatSetOption(user->Js,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); 1071 1072 ierr = MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->JsInv);CHKERRQ(ierr); 1073 ierr = MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateInvMatMult);CHKERRQ(ierr); 1074 ierr = MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateInvMatMult);CHKERRQ(ierr); 1075 ierr = MatSetOption(user->JsInv,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 1076 ierr = MatSetOption(user->JsInv,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); 1077 1078 ierr = MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 1079 ierr = MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); 1080 /* Now solve for ytrue */ 1081 ierr = KSPCreate(PETSC_COMM_WORLD,&user->solver);CHKERRQ(ierr); 1082 ierr = KSPSetFromOptions(user->solver);CHKERRQ(ierr); 1083 1084 ierr = KSPSetOperators(user->solver,user->JsBlock,user->DSG);CHKERRQ(ierr); 1085 1086 ierr = MatMult(user->JsInv,user->q,user->ytrue);CHKERRQ(ierr); 1087 /* First compute Av_u = Av*exp(-u) */ 1088 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); 1089 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); /* Note: user->u */ 1090 ierr = VecExp(user->uwork);CHKERRQ(ierr); 1091 ierr = MatMult(user->Av,user->uwork,user->Av_u);CHKERRQ(ierr); 1092 1093 /* Next update DSG = Div*S*Grad with user->u */ 1094 ierr = VecCopy(user->Av_u,user->Swork);CHKERRQ(ierr); 1095 ierr = VecReciprocal(user->Swork);CHKERRQ(ierr); 1096 if (user->use_ptap) { 1097 ierr = MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);CHKERRQ(ierr); 1098 ierr = MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);CHKERRQ(ierr); 1099 } else { 1100 ierr = MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1101 ierr = MatDiagonalScale(user->Divwork,NULL,user->Av_u);CHKERRQ(ierr); 1102 ierr = MatProductNumeric(user->DSG);CHKERRQ(ierr); 1103 } 1104 1105 /* Now solve for y */ 1106 1107 ierr = MatMult(user->JsInv,user->q,user->y);CHKERRQ(ierr); 1108 1109 user->ksp_its_initial = user->ksp_its; 1110 user->ksp_its = 0; 1111 /* Construct projection matrix Q (blocks) */ 1112 ierr = MatCreate(PETSC_COMM_WORLD,&user->Q);CHKERRQ(ierr); 1113 ierr = MatSetSizes(user->Q,dsubnlocal,ysubnlocal,user->ndata,user->ndesign);CHKERRQ(ierr); 1114 ierr = MatSetFromOptions(user->Q);CHKERRQ(ierr); 1115 ierr = MatMPIAIJSetPreallocation(user->Q,8,NULL,8,NULL);CHKERRQ(ierr); 1116 ierr = MatSeqAIJSetPreallocation(user->Q,8,NULL);CHKERRQ(ierr); 1117 1118 for (i=0; i<user->mx; i++){ 1119 x[i] = h*(i+0.5); 1120 y[i] = h*(i+0.5); 1121 z[i] = h*(i+0.5); 1122 } 1123 ierr = MatGetOwnershipRange(user->Q,&istart,&iend);CHKERRQ(ierr); 1124 1125 nx = user->mx; ny = user->mx; nz = user->mx; 1126 for (i=istart; i<iend; i++){ 1127 1128 xri = xr[i]; 1129 im = 0; 1130 xim = x[im]; 1131 while (xri>xim && im<nx){ 1132 im = im+1; 1133 xim = x[im]; 1134 } 1135 indx1 = im-1; 1136 indx2 = im; 1137 dx1 = xri - x[indx1]; 1138 dx2 = x[indx2] - xri; 1139 1140 yri = yr[i]; 1141 im = 0; 1142 yim = y[im]; 1143 while (yri>yim && im<ny){ 1144 im = im+1; 1145 yim = y[im]; 1146 } 1147 indy1 = im-1; 1148 indy2 = im; 1149 dy1 = yri - y[indy1]; 1150 dy2 = y[indy2] - yri; 1151 1152 zri = zr[i]; 1153 im = 0; 1154 zim = z[im]; 1155 while (zri>zim && im<nz){ 1156 im = im+1; 1157 zim = z[im]; 1158 } 1159 indz1 = im-1; 1160 indz2 = im; 1161 dz1 = zri - z[indz1]; 1162 dz2 = z[indz2] - zri; 1163 1164 Dx = x[indx2] - x[indx1]; 1165 Dy = y[indy2] - y[indy1]; 1166 Dz = z[indz2] - z[indz1]; 1167 1168 j = indx1 + indy1*nx + indz1*nx*ny; 1169 v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz); 1170 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1171 1172 j = indx1 + indy1*nx + indz2*nx*ny; 1173 v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz); 1174 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1175 1176 j = indx1 + indy2*nx + indz1*nx*ny; 1177 v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz); 1178 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1179 1180 j = indx1 + indy2*nx + indz2*nx*ny; 1181 v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz); 1182 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1183 1184 j = indx2 + indy1*nx + indz1*nx*ny; 1185 v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz); 1186 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1187 1188 j = indx2 + indy1*nx + indz2*nx*ny; 1189 v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz); 1190 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1191 1192 j = indx2 + indy2*nx + indz1*nx*ny; 1193 v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz); 1194 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1195 1196 j = indx2 + indy2*nx + indz2*nx*ny; 1197 v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz); 1198 ierr = MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 1199 } 1200 1201 ierr = MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1202 ierr = MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1203 /* Create MQ (composed of blocks of Q */ 1204 ierr = MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->MQ);CHKERRQ(ierr); 1205 ierr = MatShellSetOperation(user->MQ,MATOP_MULT,(void(*)(void))QMatMult);CHKERRQ(ierr); 1206 ierr = MatShellSetOperation(user->MQ,MATOP_MULT_TRANSPOSE,(void(*)(void))QMatMultTranspose);CHKERRQ(ierr); 1207 1208 /* Add noise to the measurement data */ 1209 ierr = VecSet(user->ywork,1.0);CHKERRQ(ierr); 1210 ierr = VecAYPX(user->ywork,user->noise,user->ytrue);CHKERRQ(ierr); 1211 ierr = MatMult(user->MQ,user->ywork,user->d);CHKERRQ(ierr); 1212 1213 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */ 1214 ierr = PetscFree(x);CHKERRQ(ierr); 1215 ierr = PetscFree(y);CHKERRQ(ierr); 1216 ierr = PetscFree(z);CHKERRQ(ierr); 1217 ierr = PetscLogStagePop();CHKERRQ(ierr); 1218 PetscFunctionReturn(0); 1219 } 1220 1221 PetscErrorCode EllipticDestroy(AppCtx *user) 1222 { 1223 PetscErrorCode ierr; 1224 PetscInt i; 1225 1226 PetscFunctionBegin; 1227 ierr = MatDestroy(&user->DSG);CHKERRQ(ierr); 1228 ierr = KSPDestroy(&user->solver);CHKERRQ(ierr); 1229 ierr = MatDestroy(&user->Q);CHKERRQ(ierr); 1230 ierr = MatDestroy(&user->MQ);CHKERRQ(ierr); 1231 if (!user->use_ptap) { 1232 ierr = MatDestroy(&user->Div);CHKERRQ(ierr); 1233 ierr = MatDestroy(&user->Divwork);CHKERRQ(ierr); 1234 } else { 1235 ierr = MatDestroy(&user->Diag);CHKERRQ(ierr); 1236 } 1237 if (user->use_lrc) { 1238 ierr = MatDestroy(&user->Ones);CHKERRQ(ierr); 1239 } 1240 1241 ierr = MatDestroy(&user->Grad);CHKERRQ(ierr); 1242 ierr = MatDestroy(&user->Av);CHKERRQ(ierr); 1243 ierr = MatDestroy(&user->Avwork);CHKERRQ(ierr); 1244 ierr = MatDestroy(&user->L);CHKERRQ(ierr); 1245 ierr = MatDestroy(&user->Js);CHKERRQ(ierr); 1246 ierr = MatDestroy(&user->Jd);CHKERRQ(ierr); 1247 ierr = MatDestroy(&user->JsBlock);CHKERRQ(ierr); 1248 ierr = MatDestroy(&user->JsInv);CHKERRQ(ierr); 1249 1250 ierr = VecDestroy(&user->x);CHKERRQ(ierr); 1251 ierr = VecDestroy(&user->u);CHKERRQ(ierr); 1252 ierr = VecDestroy(&user->uwork);CHKERRQ(ierr); 1253 ierr = VecDestroy(&user->utrue);CHKERRQ(ierr); 1254 ierr = VecDestroy(&user->y);CHKERRQ(ierr); 1255 ierr = VecDestroy(&user->ywork);CHKERRQ(ierr); 1256 ierr = VecDestroy(&user->ytrue);CHKERRQ(ierr); 1257 ierr = VecDestroy(&user->c);CHKERRQ(ierr); 1258 ierr = VecDestroy(&user->cwork);CHKERRQ(ierr); 1259 ierr = VecDestroy(&user->ur);CHKERRQ(ierr); 1260 ierr = VecDestroy(&user->q);CHKERRQ(ierr); 1261 ierr = VecDestroy(&user->d);CHKERRQ(ierr); 1262 ierr = VecDestroy(&user->dwork);CHKERRQ(ierr); 1263 ierr = VecDestroy(&user->lwork);CHKERRQ(ierr); 1264 ierr = VecDestroy(&user->S);CHKERRQ(ierr); 1265 ierr = VecDestroy(&user->Swork);CHKERRQ(ierr); 1266 ierr = VecDestroy(&user->Sdiag);CHKERRQ(ierr); 1267 ierr = VecDestroy(&user->Ywork);CHKERRQ(ierr); 1268 ierr = VecDestroy(&user->Twork);CHKERRQ(ierr); 1269 ierr = VecDestroy(&user->Av_u);CHKERRQ(ierr); 1270 ierr = VecDestroy(&user->js_diag);CHKERRQ(ierr); 1271 ierr = ISDestroy(&user->s_is);CHKERRQ(ierr); 1272 ierr = ISDestroy(&user->d_is);CHKERRQ(ierr); 1273 ierr = VecDestroy(&user->suby);CHKERRQ(ierr); 1274 ierr = VecDestroy(&user->subd);CHKERRQ(ierr); 1275 ierr = VecDestroy(&user->subq);CHKERRQ(ierr); 1276 ierr = VecScatterDestroy(&user->state_scatter);CHKERRQ(ierr); 1277 ierr = VecScatterDestroy(&user->design_scatter);CHKERRQ(ierr); 1278 for (i=0;i<user->ns;i++) { 1279 ierr = VecScatterDestroy(&user->yi_scatter[i]);CHKERRQ(ierr); 1280 ierr = VecScatterDestroy(&user->di_scatter[i]);CHKERRQ(ierr); 1281 } 1282 ierr = PetscFree(user->yi_scatter);CHKERRQ(ierr); 1283 ierr = PetscFree(user->di_scatter);CHKERRQ(ierr); 1284 if (user->use_lrc) { 1285 ierr = PetscFree(user->ones);CHKERRQ(ierr); 1286 ierr = MatDestroy(&user->Ones);CHKERRQ(ierr); 1287 } 1288 PetscFunctionReturn(0); 1289 } 1290 1291 PetscErrorCode EllipticMonitor(Tao tao, void *ptr) 1292 { 1293 PetscErrorCode ierr; 1294 Vec X; 1295 PetscReal unorm,ynorm; 1296 AppCtx *user = (AppCtx*)ptr; 1297 1298 PetscFunctionBegin; 1299 ierr = TaoGetSolutionVector(tao,&X);CHKERRQ(ierr); 1300 ierr = Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr); 1301 ierr = VecAXPY(user->ywork,-1.0,user->ytrue);CHKERRQ(ierr); 1302 ierr = VecAXPY(user->uwork,-1.0,user->utrue);CHKERRQ(ierr); 1303 ierr = VecNorm(user->uwork,NORM_2,&unorm);CHKERRQ(ierr); 1304 ierr = VecNorm(user->ywork,NORM_2,&ynorm);CHKERRQ(ierr); 1305 ierr = PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);CHKERRQ(ierr); 1306 PetscFunctionReturn(0); 1307 } 1308 1309 1310 1311 1312 /*TEST 1313 1314 build: 1315 requires: !complex 1316 1317 test: 1318 args: -tao_cmonitor -ns 1 -tao_type lcl -tao_gatol 1.e-3 -tao_max_it 11 1319 requires: !single 1320 1321 test: 1322 suffix: 2 1323 args: -tao_cmonitor -tao_type lcl -tao_max_it 11 -use_ptap -use_lrc -ns 1 -tao_gatol 1.e-3 1324 requires: !single 1325 1326 TEST*/ 1327