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