Actual source code: ex18.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "Isogeometric analysis of isothermal Navier-Stokes-Korteweg in 2D.";
3: #include <petscts.h>
4: #include <petscdmiga.h>
6: #define SQ(x) ((x)*(x))
8: typedef struct {
9: DM iga;
10: PetscScalar L0,h;
11: PetscScalar Ca,alpha,theta,Re;
13: // bubble centers
14: PetscScalar C1x,C1y,C2x,C2y,C3x,C3y;
15: PetscScalar R1,R2,R3;
17: } AppCtx;
19: typedef struct {
20: PetscScalar rho,ux,uy;
21: } Field;
23: PetscErrorCode InterpolateSolution(double **basis2D,Field **x,Field **xdot,PetscInt px,PetscInt py,PetscInt boffsetX,PetscInt boffsetY,
24: PetscScalar *rho,PetscScalar *rho_t,PetscScalar *rho_x,PetscScalar *rho_y,
25: PetscScalar *rho_xx,PetscScalar *rho_yy,PetscScalar *rho_xy,
26: PetscScalar *ux,PetscScalar *ux_t,PetscScalar *ux_x,PetscScalar *ux_y,
27: PetscScalar *uy,PetscScalar *uy_t,PetscScalar *uy_x,PetscScalar *uy_y);
28: PetscErrorCode FormResidual(TS ts,PetscReal t,Vec U,Vec Udot,Vec R,void *ctx);
29: PetscErrorCode FormResidualLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,Field **r,AppCtx *user);
30: PetscErrorCode FormTangent(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag,void *ctx);
31: PetscErrorCode FormTangentLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,PetscReal shift,Mat *A,AppCtx *user);
32: PetscErrorCode FormInitialCondition(AppCtx *user,Vec U);
33: PetscErrorCode WriteSolution(Vec U, const char pattern[],int number);
34: PetscErrorCode OutputMonitor(TS ts,PetscInt it_number,PetscReal c_time,Vec U,void *mctx);
38: int main(int argc, char *argv[]) {
39: PetscErrorCode ierr;
40: PetscMPIInt rank;
41: AppCtx user;
42: PetscInt p=2,N=64,C=1;
43: PetscInt ng = p+2; /* integration in each direction */
44: PetscInt Nx,Ny;
45: Vec U; /* solution vector */
46: Mat J;
47: TS ts;
48: PetscInt steps;
49: PetscReal ftime;
51: /* This code solve the dimensionless form of the isothermal
52: Navier-Stokes-Korteweg equations as presented in:
54: Gomez, Hughes, Nogueira, Calo
55: Isogeometric analysis of the isothermal Navier-Stokes-Korteweg equations
56: CMAME, 2010
58: Equation/section numbers reflect this publication.
59: */
61: // Petsc Initialization rite of passage
62: PetscInitialize(&argc,&argv,0,help);
63: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
65: // Define simulation specific parameters
66: user.L0 = 1.0; // length scale
67: user.C1x = 0.75; user.C1y = 0.50; // bubble centers
68: user.C2x = 0.25; user.C2y = 0.50;
69: user.C3x = 0.40; user.C3y = 0.75;
70: user.R1 = 0.10; user.R2 = 0.15; user.R3 = 0.08; // bubble radii
72: user.alpha = 2.0; // (Eq. 41)
73: user.theta = 0.85; // temperature parameter (just before section 5.1)
75: // Set discretization options
76: PetscOptionsBegin(PETSC_COMM_WORLD, "", "NavierStokesKorteweg Options", "IGA");
77: PetscOptionsInt("-p", "polynomial order", __FILE__, p, &p, PETSC_NULL);
78: PetscOptionsInt("-C", "global continuity order", __FILE__, C, &C, PETSC_NULL);
79: PetscOptionsInt("-N", "number of elements (along one dimension)", __FILE__, N, &N, PETSC_NULL);
80: PetscOptionsEnd();
82: // Compute simulation parameters
83: user.h = user.L0/N; // characteristic length scale of mesh (Eq. 43, simplified for uniform elements)
84: user.Ca = user.h/user.L0; // capillarity number (Eq. 38)
85: user.Re = user.alpha/user.Ca; // Reynolds number (Eq. 39)
87: // Test C < p
88: if(p <= C){
89: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Discretization inconsistent: polynomial order must be greater than degree of continuity");
90: }
92: Nx=Ny=N;
94: // Initialize B-spline space
95: DMCreate(PETSC_COMM_WORLD,&user.iga);
96: DMSetType(user.iga, DMIGA);
97: DMIGAInitializeUniform2d(user.iga,PETSC_FALSE,2,3,
98: p,Nx,C,0.0,1.0,PETSC_TRUE,ng,
99: p,Ny,C,0.0,1.0,PETSC_TRUE,ng);
101: DMCreateGlobalVector(user.iga,&U);
102: FormInitialCondition(&user,U);
103: DMIGASetFieldName(user.iga, 0, "density");
104: DMIGASetFieldName(user.iga, 1, "velocity-u");
105: DMIGASetFieldName(user.iga, 2, "velocity-v");
107: DMCreateMatrix(user.iga, MATAIJ, &J);
109: TSCreate(PETSC_COMM_WORLD,&ts);
110: TSSetType(ts,TSALPHA);
111: TSAlphaSetRadius(ts,0.5);
112: TSSetDM(ts,user.iga);
113: TSSetSolution(ts,U);
114: TSSetProblemType(ts,TS_NONLINEAR);
115: TSSetIFunction(ts,PETSC_NULL,FormResidual,&user);
116: TSSetIJacobian(ts,J,J,FormTangent,&user);
117: TSSetDuration(ts,1000000,1000.0);
118: TSSetInitialTimeStep(ts,0.0,0.001);
119: TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,PETSC_NULL);
120: TSMonitorSet(ts,OutputMonitor,PETSC_NULL,PETSC_NULL);
121: TSSetFromOptions(ts);
123: TSSolve(ts,U,&ftime);
124: TSGetTimeStepNumber(ts,&steps);
126: // Cleanup
127: TSDestroy(&ts);
128: DMDestroy(&user.iga);
129: PetscFinalize();
131: return 0;
132: }
136: PetscErrorCode FormInitialCondition(AppCtx *user,Vec U)
137: {
138: DMDALocalInfo info;
139: Field **u;
140: PetscScalar x,y,d1,d2,d3;
141: PetscInt i,j;
145: DMIGAGetLocalInfo(user->iga,&info);
146: DMIGAVecGetArray(user->iga,U,&u);
148: for(i=info.xs;i<info.xs+info.xm;i++){
149: x = user->L0*( (PetscScalar)i/(PetscScalar)info.mx );
150: for(j=info.ys;j<info.ys+info.ym;j++){
151: y = user->L0*( (PetscScalar)j/(PetscScalar)info.my );
153: d1 = PetscSqrtReal(SQ(x-user->C1x)+SQ(y-user->C1y));
154: d2 = PetscSqrtReal(SQ(x-user->C2x)+SQ(y-user->C2y));
155: d3 = PetscSqrtReal(SQ(x-user->C3x)+SQ(y-user->C3y));
157: u[j][i].rho = -0.15+0.25*( tanh(0.5*(d1-user->R1)/user->Ca) +
158: tanh(0.5*(d2-user->R2)/user->Ca) +
159: tanh(0.5*(d3-user->R3)/user->Ca) );
160: u[j][i].ux = 0.0;
161: u[j][i].uy = 0.0;
162: }
163: }
164: DMIGAVecRestoreArray(user->iga,U,&u);
165: return(0);
166: }
170: PetscErrorCode FormResidual(TS ts,PetscReal t,Vec U,Vec Udot,Vec R,void *user)
171: {
176: DMDALocalInfo info;
177: DM dm;
178: Vec localU,localUdot,localR; // local versions
179: Field **h,**hdot,**r;
181: /* get the da from the snes */
182: TSGetDM(ts, &dm);
184: /* handle the vec U */
185: DMGetLocalVector(dm,&localU);
186: DMGlobalToLocalBegin(dm,U,INSERT_VALUES,localU);
187: DMGlobalToLocalEnd(dm,U,INSERT_VALUES,localU);
189: /* handle the vec Udot */
190: DMGetLocalVector(dm,&localUdot);
191: DMGlobalToLocalBegin(dm,Udot,INSERT_VALUES,localUdot);
192: DMGlobalToLocalEnd(dm,Udot,INSERT_VALUES,localUdot);
194: /* handle the vec R */
195: DMGetLocalVector(dm,&localR);
196: VecZeroEntries(localR);
198: /* Get the arrays from the vectors */
199: DMIGAVecGetArray(dm,localU,&h);
200: DMIGAVecGetArray(dm,localUdot,&hdot);
201: DMIGAVecGetArray(dm,localR,&r);
203: /* Grab the local info and call the local residual routine */
204: DMIGAGetLocalInfo(dm,&info);
205: FormResidualLocal(&info,t,h,hdot,r,(AppCtx *) user);
207: /* Restore the arrays */
208: DMIGAVecRestoreArray(dm,localR,&r);
209: DMIGAVecRestoreArray(dm,localUdot,&hdot);
210: DMIGAVecRestoreArray(dm,localU,&h);
212: /* Add contributions back to global R */
213: VecZeroEntries(R);
214: DMLocalToGlobalBegin(dm,localR,ADD_VALUES,R); // this one adds the values
215: DMLocalToGlobalEnd(dm,localR,ADD_VALUES,R);
217: /* Restore the local vectors */
218: DMRestoreLocalVector(dm,&localU);
219: DMRestoreLocalVector(dm,&localUdot);
220: DMRestoreLocalVector(dm,&localR);
222: return(0);
223: }
227: PetscErrorCode FormResidualLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,Field **r,AppCtx *user)
228: {
229: DM iga = user->iga;
230: BD bdX, bdY;
231: PetscInt px, py, ngx, ngy;
235: DMIGAGetPolynomialOrder(iga, &px, &py, PETSC_NULL);
236: DMIGAGetNumQuadraturePoints(iga, &ngx, &ngy, PETSC_NULL);
237: DMIGAGetBasisData(iga, &bdX, &bdY, PETSC_NULL);
239: // begin and end elements for this processor
240: PetscInt bex = bdX->own_b, eex = bdX->own_e;
241: PetscInt bey = bdY->own_b, eey = bdY->own_e;
243: // Allocate space for the 3D basis to be formed
244: PetscInt Nl = (px+1)*(py+1); // number of local basis functions
245: int numD = 6; // [0] = N, [1] = dN/dx, [2] = dN/dy
246: double **basis2D;
247: ierr= PetscMalloc(numD*sizeof(double*), &basis2D);
248: int i;
249: for(i=0;i<numD;i++) {
250: PetscMalloc(Nl*sizeof(double), &basis2D[i]);
251: }
253: PetscInt ind; // temporary index variable
254: PetscInt ie,je; // iterators for elements
255: PetscInt boffsetX,boffsetY; // offsets for l2g mapping
256: PetscInt ig,jg; // iterators for gauss points
257: PetscScalar gx,gy; // gauss point locations
258: PetscScalar wgtx,wgty,wgt; // gauss point weights
259: PetscInt iba,jba; // iterators for local basis (a, matrix rows)
260: PetscScalar Nx,dNx,dNxx,Ny,dNy,dNyy; // 1D basis functions and derivatives
261: PetscInt Ax,Ay; // global matrix row/col index
262: PetscScalar Na,Na_x,Na_y,Na_xx,Na_yy,Na_xy; // 2D basis for row loop (a)
264: PetscScalar R_rho,R_ux,R_uy;
265: PetscScalar rho,rho_t,rho_x,rho_y,rho_xx,rho_yy,rho_xy;
266: PetscScalar ux,ux_t,ux_x,ux_y;
267: PetscScalar uy,uy_t,uy_x,uy_y;
268: PetscScalar tau_xx,tau_xy,tau_yx,tau_yy;
269: PetscScalar p;
271: PetscScalar Ca2 = user->Ca*user->Ca;
272: PetscScalar rRe = 1.0/user->Re;
274: for(ie=bex;ie<=eex;ie++) { // Loop over elements
275: for(je=bey;je<=eey;je++) {
277: // get basis offsets used in the local-->global mapping
278: BDGetBasisOffset(bdX,ie,&boffsetX);
279: BDGetBasisOffset(bdY,je,&boffsetY);
281: for(ig=0;ig<ngx;ig++) { // Loop over gauss points
282: for(jg=0;jg<ngy;jg++) {
284: // Get gauss point locations and weights
285: // NOTE: gauss point and weight already mapped to the parameter space
286: BDGetGaussPt(bdX,ie,ig,&gx);
287: BDGetGaussPt(bdY,je,jg,&gy);
288: BDGetGaussWt(bdX,ie,ig,&wgtx);
289: BDGetGaussWt(bdY,je,jg,&wgty);
291: wgt = wgtx*wgty;
293: for(jba=0;jba<(py+1);jba++) { // Assemble the 2D basis
294: for(iba=0;iba<(px+1);iba++) {
296: BDGetBasis(bdX,ie,ig,iba,0,&Nx);
297: BDGetBasis(bdX,ie,ig,iba,1,&dNx);
298: BDGetBasis(bdX,ie,ig,iba,2,&dNxx);
300: BDGetBasis(bdY,je,jg,jba,0,&Ny);
301: BDGetBasis(bdY,je,jg,jba,1,&dNy);
302: BDGetBasis(bdY,je,jg,jba,2,&dNyy);
304: // 2D basis is a tensor product
305: ind = jba*(px+1)+iba;
306: basis2D[0][ind] = Nx * Ny; // N
307: basis2D[1][ind] = dNx * Ny; // dN/dx
308: basis2D[2][ind] = Nx * dNy; // dN/dy
309: basis2D[3][ind] = dNxx * Ny; // d^2N/dx^2
310: basis2D[4][ind] = Nx * dNyy; // d^2N/dy^2
311: basis2D[5][ind] = dNx * dNy; // d^2N/dxdy
313: }
314: } // end 2D basis assembly
316: // Problem coefficient evaluation
317: InterpolateSolution(basis2D,h,hdot,px,py,boffsetX,boffsetY,
318: &rho,&rho_t,&rho_x,&rho_y,
319: &rho_xx,&rho_yy,&rho_xy,
320: &ux,&ux_t,&ux_x,&ux_y,
321: &uy,&uy_t,&uy_x,&uy_y);
323: // compute pressure (Eq. 34.3)
324: p = 8.0/27.0*user->theta*rho/(1.0-rho)-rho*rho;
326: // compute viscous stress tensor (Eq. 34.4)
327: tau_xx = 2.0*ux_x - 2.0/3.0*(ux_x+uy_y);
328: tau_xy = ux_y + uy_x ;
329: tau_yy = 2.0*uy_y - 2.0/3.0*(ux_x+uy_y);
330: tau_yx = tau_xy;
332: for(jba=0;jba<(py+1);jba++) { // loop over basis 1st time (a, matrix rows)
333: for(iba=0;iba<(px+1);iba++) {
335: Ax = boffsetX+iba; // local to global map
336: Ay = boffsetY+jba;
338: ind = jba*(px+1)+iba;
339: Na = basis2D[0][ind];
340: Na_x = basis2D[1][ind];
341: Na_y = basis2D[2][ind];
342: Na_xx = basis2D[3][ind];
343: Na_yy = basis2D[4][ind];
344: Na_xy = basis2D[5][ind];
346: // (Eq. 19, modified to be dimensionless)
347: R_rho = Na*rho_t;
348: R_rho += -rho*(Na_x*ux + Na_y*uy);
350: R_ux = Na*ux*rho_t;
351: R_ux += Na*rho*ux_t;
352: R_ux += -rho*(Na_x*ux*ux + Na_y*ux*uy);
353: R_ux += -Na_x*p;
354: R_ux += rRe*(Na_x*tau_xx + Na_y*tau_xy);
355: //R_ux += -Ca2*rho*rho_x*(Na_xx+Na_xy);
356: //R_ux += -Ca2*Na_x*(rho_x*rho_x+rho_y*rho_y);
357: //R_ux += -Ca2*(rho_xx*Na + rho_x*Na_x + rho_xy*Na + rho_y*Na_x)*rho_x;
358: // rewritten uses Victor's corrections
359: R_ux += -Ca2*(Na_xx*rho_x + Na_xy*rho_y);
360: R_ux += -Ca2*Na_x*(rho_x*rho_x+rho_y*rho_y);
361: R_ux += -Ca2*Na*(rho_xx*rho_x+rho_xy*rho_y);
362: R_ux += -Ca2*rho_x*(Na_x*rho_x+Na_y*rho_y);
364: R_uy = Na*uy*rho_t;
365: R_uy += Na*rho*uy_t;
366: R_uy += -rho*(Na_x*uy*ux + Na_y*uy*uy);
367: R_uy += -Na_y*p;
368: R_uy += rRe*(Na_x*tau_yx + Na_y*tau_yy);
370: //R_uy += -Ca2*rho*rho_y*(Na_xy+Na_yy);
371: //R_uy += -Ca2*Na_y*(rho_x*rho_x+rho_y*rho_y);
372: //R_uy += -Ca2*(rho_xy*Na + rho_x*Na_y + rho_yy*Na + rho_y*Na_y)*rho_y;
374: R_uy += -Ca2*(Na_xy*rho_x + Na_yy*rho_y);
375: R_uy += -Ca2*Na_y*(rho_x*rho_x+rho_y*rho_y);
376: R_uy += -Ca2*Na*(rho_xy*rho_x+rho_yy*rho_y);
377: R_uy += -Ca2*rho_y*(Na_x*rho_x+Na_y*rho_y);
379: r[Ay][Ax].rho += R_rho*wgt;
380: r[Ay][Ax].ux += R_ux*wgt;
381: r[Ay][Ax].uy += R_uy*wgt;
383: }
384: } // end basis a loop
386: }
387: } // end gauss point loop
389: }
390: } // end element loop
392: for(i=0;i<numD;i++) {
393: PetscFree(basis2D[i]);
394: }
395: PetscFree(basis2D);
398: return(0);
399: }
403: PetscErrorCode FormTangent(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag,void *ctx)
404: {
408: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "FormTangent not implemented, use -snes_mf");
410: DMDALocalInfo info;
411: DM da_dof;
412: Vec localU,localUdot; // local versions
413: Field **h,**hdot;
415: /* get the da from the snes */
416: TSGetDM(ts,(DM*)&da_dof);
418: /* handle the vec U */
419: DMGetLocalVector(da_dof,&localU);
420: DMGlobalToLocalBegin(da_dof,U,INSERT_VALUES,localU);
421: DMGlobalToLocalEnd(da_dof,U,INSERT_VALUES,localU);
423: /* handle the vec Udot */
424: DMGetLocalVector(da_dof,&localUdot);
425: DMGlobalToLocalBegin(da_dof,Udot,INSERT_VALUES,localUdot);
426: DMGlobalToLocalEnd(da_dof,Udot,INSERT_VALUES,localUdot);
428: /* Get the arrays from the vectors */
429: DMDAVecGetArray(da_dof,localU,&h);
430: DMDAVecGetArray(da_dof,localUdot,&hdot);
432: /* Grab the local info and call the local tangent routine */
433: DMDAGetLocalInfo(da_dof,&info);
434: MatZeroEntries(*B); // pre-zero the matrix
435: FormTangentLocal(&info,t,h,hdot,shift,B,(AppCtx *) ctx);
436: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
437: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
438: if (*A != *B) { // then we could be matrix free
439: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
440: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
441: }
442: *flag = SAME_NONZERO_PATTERN; /* the sparsity pattern does not change */
444: DMDAVecRestoreArray(da_dof,localUdot,&hdot);
445: DMDAVecRestoreArray(da_dof,localU,&h);
447: /* Restore the arrays and local vectors */
448: DMRestoreLocalVector(da_dof,&localU);
449: DMRestoreLocalVector(da_dof,&localUdot);
451: return(0);
452: }
456: PetscErrorCode FormTangentLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,PetscReal shift,Mat *A,AppCtx *user)
457: {
460: return(0);
461: }
465: PetscErrorCode InterpolateSolution(double **basis2D,Field **x,Field **xdot,PetscInt px,PetscInt py,PetscInt boffsetX,PetscInt boffsetY,
466: PetscScalar *rho,PetscScalar *rho_t,PetscScalar *rho_x,PetscScalar *rho_y,
467: PetscScalar *rho_xx,PetscScalar *rho_yy,PetscScalar *rho_xy,
468: PetscScalar *ux,PetscScalar *ux_t,PetscScalar *ux_x,PetscScalar *ux_y,
469: PetscScalar *uy,PetscScalar *uy_t,PetscScalar *uy_x,PetscScalar *uy_y)
470: {
472: (*rho) = 0.0; (*rho_x) = 0.0; (*rho_y) = 0.0; (*rho_t) = 0.0;
473: (*rho_xx) = 0.0; (*rho_yy) = 0.0; (*rho_xy) = 0.0;
474: (*ux) = 0.0; (*ux_x) = 0.0; (*ux_y) = 0.0; (*ux_t) = 0.0;
475: (*uy) = 0.0; (*uy_x) = 0.0; (*uy_y) = 0.0; (*uy_t) = 0.0;
477: int ipa,jpa,ind;
478: for(jpa=0;jpa<(py+1);jpa++) {
479: for(ipa=0;ipa<(px+1);ipa++) {
481: ind = jpa*(px+1)+ipa;
482: (*rho) += basis2D[0][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
483: (*ux) += basis2D[0][ind] * x[boffsetY+jpa][boffsetX+ipa].ux;
484: (*uy) += basis2D[0][ind] * x[boffsetY+jpa][boffsetX+ipa].uy;
486: (*rho_x) += basis2D[1][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
487: (*ux_x) += basis2D[1][ind] * x[boffsetY+jpa][boffsetX+ipa].ux;
488: (*uy_x) += basis2D[1][ind] * x[boffsetY+jpa][boffsetX+ipa].uy;
490: (*rho_y) += basis2D[2][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
491: (*ux_y) += basis2D[2][ind] * x[boffsetY+jpa][boffsetX+ipa].ux;
492: (*uy_y) += basis2D[2][ind] * x[boffsetY+jpa][boffsetX+ipa].uy;
494: (*rho_xx) += basis2D[3][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
495: (*rho_yy) += basis2D[4][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
496: (*rho_xy) += basis2D[5][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
498: (*rho_t) += basis2D[0][ind] * xdot[boffsetY+jpa][boffsetX+ipa].rho;
499: (*ux_t) += basis2D[0][ind] * xdot[boffsetY+jpa][boffsetX+ipa].ux;
500: (*uy_t) += basis2D[0][ind] * xdot[boffsetY+jpa][boffsetX+ipa].uy;
502: }
503: }
506: return(0);
507: }
512: PetscErrorCode WriteSolution(Vec U, const char pattern[],int number)
513: {
515: PetscErrorCode ierr;
516: MPI_Comm comm;
517: char filename[256];
518: PetscViewer viewer;
521: sprintf(filename,pattern,number);
522: PetscObjectGetComm((PetscObject)U,&comm);
523: PetscViewerBinaryOpen(comm,filename,FILE_MODE_WRITE,&viewer);
524: VecView(U,viewer);
525: PetscViewerFlush(viewer);
526: PetscViewerDestroy(&viewer);
527: return(0);
528: }
532: PetscErrorCode OutputMonitor(TS ts,PetscInt it_number,PetscReal c_time,Vec U,void *mctx)
533: {
537: WriteSolution(U,"./nsk%d.dat",it_number);
539: return(0);
540: }