Actual source code: ex3.c
petsc-3.8.4 2018-03-24
2: static char help[] = "Solves 1D heat equation with FEM formulation.\n\
3: Input arguments are\n\
4: -useAlhs: solve Alhs*U' = (Arhs*U + g) \n\
5: otherwise, solve U' = inv(Alhs)*(Arhs*U + g) \n\n";
7: /*--------------------------------------------------------------------------
8: Solves 1D heat equation U_t = U_xx with FEM formulation:
9: Alhs*U' = rhs (= Arhs*U + g)
10: We thank Chris Cox <clcox@clemson.edu> for contributing the original code
11: ----------------------------------------------------------------------------*/
13: #include <petscksp.h>
14: #include <petscts.h>
16: /* special variable - max size of all arrays */
17: #define num_z 10
19: /*
20: User-defined application context - contains data needed by the
21: application-provided call-back routines.
22: */
23: typedef struct {
24: Mat Amat; /* left hand side matrix */
25: Vec ksp_rhs,ksp_sol; /* working vectors for formulating inv(Alhs)*(Arhs*U+g) */
26: int max_probsz; /* max size of the problem */
27: PetscBool useAlhs; /* flag (1 indicates solving Alhs*U' = Arhs*U+g */
28: int nz; /* total number of grid points */
29: PetscInt m; /* total number of interio grid points */
30: Vec solution; /* global exact ts solution vector */
31: PetscScalar *z; /* array of grid points */
32: PetscBool debug; /* flag (1 indicates activation of debugging printouts) */
33: } AppCtx;
35: extern PetscScalar exact(PetscScalar,PetscReal);
36: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
37: extern PetscErrorCode Petsc_KSPSolve(AppCtx*);
38: extern PetscScalar bspl(PetscScalar*,PetscScalar,PetscInt,PetscInt,PetscInt[][2],PetscInt);
39: extern PetscErrorCode femBg(PetscScalar[][3],PetscScalar*,PetscInt,PetscScalar*,PetscReal);
40: extern PetscErrorCode femA(AppCtx*,PetscInt,PetscScalar*);
41: extern PetscErrorCode rhs(AppCtx*,PetscScalar*, PetscInt, PetscScalar*,PetscReal);
42: extern PetscErrorCode RHSfunction(TS,PetscReal,Vec,Vec,void*);
44: int main(int argc,char **argv)
45: {
46: PetscInt i,m,nz,steps,max_steps,k,nphase=1;
47: PetscScalar zInitial,zFinal,val,*z;
48: PetscReal stepsz[4],T,ftime;
50: TS ts;
51: SNES snes;
52: Mat Jmat;
53: AppCtx appctx; /* user-defined application context */
54: Vec init_sol; /* ts solution vector */
55: PetscMPIInt size;
57: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
58: MPI_Comm_size(PETSC_COMM_WORLD,&size);
59: if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only");
61: PetscOptionsGetInt(NULL,NULL,"-nphase",&nphase,NULL);
62: if (nphase > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nphase must be an integer between 1 and 3");
64: /* initializations */
65: zInitial = 0.0;
66: zFinal = 1.0;
67: T = 0.014/nphase;
68: nz = num_z;
69: m = nz-2;
70: appctx.nz = nz;
71: max_steps = (PetscInt)10000;
73: appctx.m = m;
74: appctx.max_probsz = nz;
75: appctx.debug = PETSC_FALSE;
76: appctx.useAlhs = PETSC_FALSE;
78: PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug);
79: PetscOptionsHasName(NULL,NULL,"-useAlhs",&appctx.useAlhs);
81: /* create vector to hold ts solution */
82: /*-----------------------------------*/
83: VecCreate(PETSC_COMM_WORLD, &init_sol);
84: VecSetSizes(init_sol, PETSC_DECIDE, m);
85: VecSetFromOptions(init_sol);
87: /* create vector to hold true ts soln for comparison */
88: VecDuplicate(init_sol, &appctx.solution);
90: /* create LHS matrix Amat */
91: /*------------------------*/
92: MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat);
93: MatSetFromOptions(appctx.Amat);
94: MatSetUp(appctx.Amat);
95: /* set space grid points - interio points only! */
96: PetscMalloc1(nz+1,&z);
97: for (i=0; i<nz; i++) z[i]=(i)*((zFinal-zInitial)/(nz-1));
98: appctx.z = z;
99: femA(&appctx,nz,z);
101: /* create the jacobian matrix */
102: /*----------------------------*/
103: MatCreate(PETSC_COMM_WORLD, &Jmat);
104: MatSetSizes(Jmat,PETSC_DECIDE,PETSC_DECIDE,m,m);
105: MatSetFromOptions(Jmat);
106: MatSetUp(Jmat);
108: /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
109: VecDuplicate(init_sol,&appctx.ksp_rhs);
110: VecDuplicate(init_sol,&appctx.ksp_sol);
112: /* set intial guess */
113: /*------------------*/
114: for (i=0; i<nz-2; i++) {
115: val = exact(z[i+1], 0.0);
116: VecSetValue(init_sol,i,(PetscScalar)val,INSERT_VALUES);
117: }
118: VecAssemblyBegin(init_sol);
119: VecAssemblyEnd(init_sol);
121: /*create a time-stepping context and set the problem type */
122: /*--------------------------------------------------------*/
123: TSCreate(PETSC_COMM_WORLD, &ts);
124: TSSetProblemType(ts,TS_NONLINEAR);
126: /* set time-step method */
127: TSSetType(ts,TSCN);
129: /* Set optional user-defined monitoring routine */
130: TSMonitorSet(ts,Monitor,&appctx,NULL);
131: /* set the right hand side of U_t = RHSfunction(U,t) */
132: TSSetRHSFunction(ts,NULL,(PetscErrorCode (*)(TS,PetscScalar,Vec,Vec,void*))RHSfunction,&appctx);
134: if (appctx.useAlhs) {
135: /* set the left hand side matrix of Amat*U_t = rhs(U,t) */
137: /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the
138: * Alhs matrix without making a copy. Either finite difference the entire thing or use analytic Jacobians in both
139: * places.
140: */
141: TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,&appctx);
142: TSSetIJacobian(ts,appctx.Amat,appctx.Amat,TSComputeIJacobianConstant,&appctx);
143: }
145: /* use petsc to compute the jacobian by finite differences */
146: TSGetSNES(ts,&snes);
147: SNESSetJacobian(snes,Jmat,Jmat,SNESComputeJacobianDefault,NULL);
149: /* get the command line options if there are any and set them */
150: TSSetFromOptions(ts);
152: #if defined(PETSC_HAVE_SUNDIALS)
153: {
154: TSType type;
155: PetscBool sundialstype=PETSC_FALSE;
156: TSGetType(ts,&type);
157: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundialstype);
158: if (sundialstype && appctx.useAlhs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use Alhs formulation for TSSUNDIALS type");
159: }
160: #endif
161: /* Sets the initial solution */
162: TSSetSolution(ts,init_sol);
164: stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */
165: ftime = 0.0;
166: for (k=0; k<nphase; k++) {
167: if (nphase > 1) {PetscPrintf(PETSC_COMM_WORLD,"Phase %D initial time %g, stepsz %g, duration: %g\n",k,(double)ftime,(double)stepsz[k],(double)((k+1)*T));}
168: TSSetTime(ts,ftime);
169: TSSetTimeStep(ts,stepsz[k]);
170: TSSetMaxSteps(ts,max_steps);
171: TSSetMaxTime(ts,(k+1)*T);
172: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
174: /* loop over time steps */
175: /*----------------------*/
176: TSSolve(ts,init_sol);
177: TSGetSolveTime(ts,&ftime);
178: TSGetStepNumber(ts,&steps);
179: stepsz[k+1] = stepsz[k]*1.5; /* change step size for the next phase */
180: }
182: /* free space */
183: TSDestroy(&ts);
184: MatDestroy(&appctx.Amat);
185: MatDestroy(&Jmat);
186: VecDestroy(&appctx.ksp_rhs);
187: VecDestroy(&appctx.ksp_sol);
188: VecDestroy(&init_sol);
189: VecDestroy(&appctx.solution);
190: PetscFree(z);
192: PetscFinalize();
193: return ierr;
194: }
196: /*------------------------------------------------------------------------
197: Set exact solution
198: u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t)
199: --------------------------------------------------------------------------*/
200: PetscScalar exact(PetscScalar z,PetscReal t)
201: {
202: PetscScalar val, ex1, ex2;
204: ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t);
205: ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t);
206: val = PetscSinScalar(6*PETSC_PI*z)*ex1 + 3.*PetscSinScalar(2*PETSC_PI*z)*ex2;
207: return val;
208: }
210: /*
211: Monitor - User-provided routine to monitor the solution computed at
212: each timestep. This example plots the solution and computes the
213: error in two different norms.
215: Input Parameters:
216: ts - the timestep context
217: step - the count of the current step (with 0 meaning the
218: initial condition)
219: time - the current time
220: u - the solution at this timestep
221: ctx - the user-provided context for this monitoring routine.
222: In this case we use the application context which contains
223: information about the problem size, workspace and the exact
224: solution.
225: */
226: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
227: {
228: AppCtx *appctx = (AppCtx*)ctx;
230: PetscInt i,m=appctx->m;
231: PetscReal norm_2,norm_max,h=1.0/(m+1);
232: PetscScalar *u_exact;
234: /* Compute the exact solution */
235: VecGetArray(appctx->solution,&u_exact);
236: for (i=0; i<m; i++) u_exact[i] = exact(appctx->z[i+1],time);
237: VecRestoreArray(appctx->solution,&u_exact);
239: /* Print debugging information if desired */
240: if (appctx->debug) {
241: PetscPrintf(PETSC_COMM_SELF,"Computed solution vector at time %g\n",(double)time);
242: VecView(u,PETSC_VIEWER_STDOUT_SELF);
243: PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");
244: VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
245: }
247: /* Compute the 2-norm and max-norm of the error */
248: VecAXPY(appctx->solution,-1.0,u);
249: VecNorm(appctx->solution,NORM_2,&norm_2);
251: norm_2 = PetscSqrtReal(h)*norm_2;
252: VecNorm(appctx->solution,NORM_MAX,&norm_max);
254: PetscPrintf(PETSC_COMM_SELF,"Timestep %D: time = %g, 2-norm error = %6.4f, max norm error = %6.4f\n",step,(double)time,(double)norm_2,(double)norm_max);
256: /*
257: Print debugging information if desired
258: */
259: if (appctx->debug) {
260: PetscPrintf(PETSC_COMM_SELF,"Error vector\n");
261: VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
262: }
263: return 0;
264: }
266: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267: %% Function to solve a linear system using KSP %%
268: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
270: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
271: {
273: KSP ksp;
274: PC pc;
276: /*create the ksp context and set the operators,that is, associate the system matrix with it*/
277: KSPCreate(PETSC_COMM_WORLD,&ksp);
278: KSPSetOperators(ksp,obj->Amat,obj->Amat);
280: /*get the preconditioner context, set its type and the tolerances*/
281: KSPGetPC(ksp,&pc);
282: PCSetType(pc,PCLU);
283: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
285: /*get the command line options if there are any and set them*/
286: KSPSetFromOptions(ksp);
288: /*get the linear system (ksp) solve*/
289: KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol);
291: KSPDestroy(&ksp);
292: return 0;
293: }
295: /***********************************************************************
296: * Function to return value of basis function or derivative of basis *
297: * function. *
298: ***********************************************************************
299: * *
300: * Arguments: *
301: * x = array of xpoints or nodal values *
302: * xx = point at which the basis function is to be *
303: * evaluated. *
304: * il = interval containing xx. *
305: * iq = indicates which of the two basis functions in *
306: * interval intrvl should be used *
307: * nll = array containing the endpoints of each interval. *
308: * id = If id ~= 2, the value of the basis function *
309: * is calculated; if id = 2, the value of the *
310: * derivative of the basis function is returned. *
311: ***********************************************************************/
313: PetscScalar bspl(PetscScalar *x, PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id)
314: {
315: PetscScalar x1,x2,bfcn;
316: PetscInt i1,i2,iq1,iq2;
318: /*** Determine which basis function in interval intrvl is to be used in ***/
319: iq1 = iq;
320: if (iq1==0) iq2 = 1;
321: else iq2 = 0;
323: /*** Determine endpoint of the interval intrvl ***/
324: i1=nll[il][iq1];
325: i2=nll[il][iq2];
327: /*** Determine nodal values at the endpoints of the interval intrvl ***/
328: x1=x[i1];
329: x2=x[i2];
330: /* printf("x1=%g\tx2=%g\txx=%g\n",x1,x2,xx); */
331: /*** Evaluate basis function ***/
332: if (id == 2) bfcn=(1.0)/(x1-x2);
333: else bfcn=(xx-x2)/(x1-x2);
334: /* printf("bfcn=%g\n",bfcn); */
335: return bfcn;
336: }
338: /*---------------------------------------------------------
339: Function called by rhs function to get B and g
340: ---------------------------------------------------------*/
341: PetscErrorCode femBg(PetscScalar btri[][3],PetscScalar *f,PetscInt nz,PetscScalar *z, PetscReal t)
342: {
343: PetscInt i,j,jj,il,ip,ipp,ipq,iq,iquad,iqq;
344: PetscInt nli[num_z][2],indx[num_z];
345: PetscScalar dd,dl,zip,zipq,zz,b_z,bb_z,bij;
346: PetscScalar zquad[num_z][3],dlen[num_z],qdwt[3];
348: /* initializing everything - btri and f are initialized in rhs.c */
349: for (i=0; i < nz; i++) {
350: nli[i][0] = 0;
351: nli[i][1] = 0;
352: indx[i] = 0;
353: zquad[i][0] = 0.0;
354: zquad[i][1] = 0.0;
355: zquad[i][2] = 0.0;
356: dlen[i] = 0.0;
357: } /*end for (i)*/
359: /* quadrature weights */
360: qdwt[0] = 1.0/6.0;
361: qdwt[1] = 4.0/6.0;
362: qdwt[2] = 1.0/6.0;
364: /* 1st and last nodes have Dirichlet boundary condition -
365: set indices there to -1 */
367: for (i=0; i < nz-1; i++) indx[i] = i-1;
368: indx[nz-1] = -1;
370: ipq = 0;
371: for (il=0; il < nz-1; il++) {
372: ip = ipq;
373: ipq = ip+1;
374: zip = z[ip];
375: zipq = z[ipq];
376: dl = zipq-zip;
377: zquad[il][0] = zip;
378: zquad[il][1] = (0.5)*(zip+zipq);
379: zquad[il][2] = zipq;
380: dlen[il] = PetscAbsScalar(dl);
381: nli[il][0] = ip;
382: nli[il][1] = ipq;
383: }
385: for (il=0; il < nz-1; il++) {
386: for (iquad=0; iquad < 3; iquad++) {
387: dd = (dlen[il])*(qdwt[iquad]);
388: zz = zquad[il][iquad];
390: for (iq=0; iq < 2; iq++) {
391: ip = nli[il][iq];
392: b_z = bspl(z,zz,il,iq,nli,2);
393: i = indx[ip];
395: if (i > -1) {
396: for (iqq=0; iqq < 2; iqq++) {
397: ipp = nli[il][iqq];
398: bb_z = bspl(z,zz,il,iqq,nli,2);
399: j = indx[ipp];
400: bij = -b_z*bb_z;
402: if (j > -1) {
403: jj = 1+j-i;
404: btri[i][jj] += bij*dd;
405: } else {
406: f[i] += bij*dd*exact(z[ipp], t);
407: /* f[i] += 0.0; */
408: /* if (il==0 && j==-1) { */
409: /* f[i] += bij*dd*exact(zz,t); */
410: /* }*/ /*end if*/
411: } /*end else*/
412: } /*end for (iqq)*/
413: } /*end if (i>0)*/
414: } /*end for (iq)*/
415: } /*end for (iquad)*/
416: } /*end for (il)*/
417: return 0;
418: }
420: PetscErrorCode femA(AppCtx *obj,PetscInt nz,PetscScalar *z)
421: {
422: PetscInt i,j,il,ip,ipp,ipq,iq,iquad,iqq;
423: PetscInt nli[num_z][2],indx[num_z];
424: PetscScalar dd,dl,zip,zipq,zz,bb,bbb,aij;
425: PetscScalar rquad[num_z][3],dlen[num_z],qdwt[3],add_term;
428: /* initializing everything */
430: for (i=0; i < nz; i++) {
431: nli[i][0] = 0;
432: nli[i][1] = 0;
433: indx[i] = 0;
434: rquad[i][0] = 0.0;
435: rquad[i][1] = 0.0;
436: rquad[i][2] = 0.0;
437: dlen[i] = 0.0;
438: } /*end for (i)*/
440: /* quadrature weights */
441: qdwt[0] = 1.0/6.0;
442: qdwt[1] = 4.0/6.0;
443: qdwt[2] = 1.0/6.0;
445: /* 1st and last nodes have Dirichlet boundary condition -
446: set indices there to -1 */
448: for (i=0; i < nz-1; i++) indx[i]=i-1;
449: indx[nz-1]=-1;
451: ipq = 0;
453: for (il=0; il < nz-1; il++) {
454: ip = ipq;
455: ipq = ip+1;
456: zip = z[ip];
457: zipq = z[ipq];
458: dl = zipq-zip;
459: rquad[il][0] = zip;
460: rquad[il][1] = (0.5)*(zip+zipq);
461: rquad[il][2] = zipq;
462: dlen[il] = PetscAbsScalar(dl);
463: nli[il][0] = ip;
464: nli[il][1] = ipq;
465: } /*end for (il)*/
467: for (il=0; il < nz-1; il++) {
468: for (iquad=0; iquad < 3; iquad++) {
469: dd = (dlen[il])*(qdwt[iquad]);
470: zz = rquad[il][iquad];
472: for (iq=0; iq < 2; iq++) {
473: ip = nli[il][iq];
474: bb = bspl(z,zz,il,iq,nli,1);
475: i = indx[ip];
476: if (i > -1) {
477: for (iqq=0; iqq < 2; iqq++) {
478: ipp = nli[il][iqq];
479: bbb = bspl(z,zz,il,iqq,nli,1);
480: j = indx[ipp];
481: aij = bb*bbb;
482: if (j > -1) {
483: add_term = aij*dd;
484: MatSetValue(obj->Amat,i,j,add_term,ADD_VALUES);
485: }/*endif*/
486: } /*end for (iqq)*/
487: } /*end if (i>0)*/
488: } /*end for (iq)*/
489: } /*end for (iquad)*/
490: } /*end for (il)*/
491: MatAssemblyBegin(obj->Amat,MAT_FINAL_ASSEMBLY);
492: MatAssemblyEnd(obj->Amat,MAT_FINAL_ASSEMBLY);
493: return 0;
494: }
496: /*---------------------------------------------------------
497: Function to fill the rhs vector with
498: By + g values ****
499: ---------------------------------------------------------*/
500: PetscErrorCode rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
501: {
502: PetscInt i,j,js,je,jj;
503: PetscScalar val,g[num_z],btri[num_z][3],add_term;
506: for (i=0; i < nz-2; i++) {
507: for (j=0; j <= 2; j++) btri[i][j]=0.0;
508: g[i] = 0.0;
509: }
511: /* call femBg to set the tri-diagonal b matrix and vector g */
512: femBg(btri,g,nz,z,t);
514: /* setting the entries of the right hand side vector */
515: for (i=0; i < nz-2; i++) {
516: val = 0.0;
517: js = 0;
518: if (i == 0) js = 1;
519: je = 2;
520: if (i == nz-2) je = 1;
522: for (jj=js; jj <= je; jj++) {
523: j = i+jj-1;
524: val += (btri[i][jj])*(y[j]);
525: }
526: add_term = val + g[i];
527: VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);
528: }
529: VecAssemblyBegin(obj->ksp_rhs);
530: VecAssemblyEnd(obj->ksp_rhs);
532: /* return to main driver function */
533: return 0;
534: }
536: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
537: %% Function to form the right hand side of the time-stepping problem. %%
538: %% -------------------------------------------------------------------------------------------%%
539: if (useAlhs):
540: globalout = By+g
541: else if (!useAlhs):
542: globalout = f(y,t)=Ainv(By+g),
543: in which the ksp solver to transform the problem A*ydot=By+g
544: to the problem ydot=f(y,t)=inv(A)*(By+g)
545: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
547: PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
548: {
549: PetscErrorCode ierr;
550: AppCtx *obj = (AppCtx*)ctx;
551: PetscScalar soln[num_z];
552: const PetscScalar *soln_ptr;
553: PetscInt i,nz=obj->nz;
554: PetscReal time;
556: /* get the previous solution to compute updated system */
557: VecGetArrayRead(globalin,&soln_ptr);
558: for (i=0; i < num_z-2; i++) soln[i] = soln_ptr[i];
559: VecRestoreArrayRead(globalin,&soln_ptr);
560: soln[num_z-1] = 0.0;
561: soln[num_z-2] = 0.0;
563: /* clear out the matrix and rhs for ksp to keep things straight */
564: VecSet(obj->ksp_rhs,(PetscScalar)0.0);
566: time = t;
567: /* get the updated system */
568: rhs(obj,soln,nz,obj->z,time); /* setup of the By+g rhs */
570: /* do a ksp solve to get the rhs for the ts problem */
571: if (obj->useAlhs) {
572: /* ksp_sol = ksp_rhs */
573: VecCopy(obj->ksp_rhs,globalout);
574: } else {
575: /* ksp_sol = inv(Amat)*ksp_rhs */
576: Petsc_KSPSolve(obj);
577: VecCopy(obj->ksp_sol,globalout);
578: }
579: return 0;
580: }