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