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