Actual source code: ex18.c
1: #include <petsc.h>
2: #include <petscdmiga.h>
4: #define SQ(x) ((x)*(x))
6: typedef struct {
7: DM iga;
8: PetscScalar L0,h;
9: PetscScalar Ca,alpha,theta,Re;
11: // bubble centers
12: PetscScalar C1x,C1y,C2x,C2y,C3x,C3y;
13: PetscScalar R1,R2,R3;
15: } AppCtx;
17: typedef struct {
18: PetscScalar rho,ux,uy;
19: } Field;
21: PetscErrorCode InterpolateSolution(double **basis2D,Field **x,Field **xdot,PetscInt px,PetscInt py,PetscInt boffsetX,PetscInt boffsetY,
22: PetscScalar *rho,PetscScalar *rho_t,PetscScalar *rho_x,PetscScalar *rho_y,
23: PetscScalar *rho_xx,PetscScalar *rho_yy,PetscScalar *rho_xy,
24: PetscScalar *ux,PetscScalar *ux_t,PetscScalar *ux_x,PetscScalar *ux_y,
25: PetscScalar *uy,PetscScalar *uy_t,PetscScalar *uy_x,PetscScalar *uy_y);
26: PetscErrorCode FormResidual(TS ts,PetscReal t,Vec U,Vec Udot,Vec R,void *ctx);
27: PetscErrorCode FormResidualLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,Field **r,AppCtx *user);
28: PetscErrorCode FormTangent(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag,void *ctx);
29: PetscErrorCode FormTangentLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,PetscReal shift,Mat *A,AppCtx *user);
30: PetscErrorCode FormInitialCondition(AppCtx *user,Vec U);
31: PetscErrorCode WriteSolution(Vec U, const char pattern[],int number);
32: PetscErrorCode OutputMonitor(TS ts,PetscInt it_number,PetscReal c_time,Vec U,void *mctx);
36: int main(int argc, char *argv[]) {
37: PetscErrorCode ierr;
38: PetscMPIInt rank;
39: AppCtx user;
40: PetscInt p=2,N=64,C=1;
41: PetscInt ng = p+2; /* integration in each direction */
42: PetscInt Nx,Ny;
43: Vec U; /* solution vector */
44: Mat J;
45: TS ts;
46: PetscInt steps;
47: PetscReal ftime;
49: /* This code solve the dimensionless form of the isothermal
50: Navier-Stokes-Korteweg equations as presented in:
52: Gomez, Hughes, Nogueira, Calo
53: Isogeometric analysis of the isothermal Navier-Stokes-Korteweg equations
54: CMAME, 2010
56: Equation/section numbers reflect this publication.
57: */
59: // Petsc Initialization rite of passage
60: PetscInitialize(&argc,&argv,0,0);
61: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
63: // Define simulation specific parameters
64: user.L0 = 1.0; // length scale
65: user.C1x = 0.75; user.C1y = 0.50; // bubble centers
66: user.C2x = 0.25; user.C2y = 0.50;
67: user.C3x = 0.40; user.C3y = 0.75;
68: user.R1 = 0.10; user.R2 = 0.15; user.R3 = 0.08; // bubble radii
70: user.alpha = 2.0; // (Eq. 41)
71: user.theta = 0.85; // temperature parameter (just before section 5.1)
73: // Set discretization options
74: PetscOptionsBegin(PETSC_COMM_WORLD, "", "NavierStokesKorteweg Options", "IGA");
75: PetscOptionsInt("-p", "polynomial order", __FILE__, p, &p, PETSC_NULL);
76: PetscOptionsInt("-C", "global continuity order", __FILE__, C, &C, PETSC_NULL);
77: PetscOptionsInt("-N", "number of elements (along one dimension)", __FILE__, N, &N, PETSC_NULL);
78: PetscOptionsEnd();
80: // Compute simulation parameters
81: user.h = user.L0/N; // characteristic length scale of mesh (Eq. 43, simplified for uniform elements)
82: user.Ca = user.h/user.L0; // capillarity number (Eq. 38)
83: user.Re = user.alpha/user.Ca; // Reynolds number (Eq. 39)
85: // Test C < p
86: if(p <= C){
87: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Discretization inconsistent: polynomial order must be greater than degree of continuity");
88: }
90: Nx=Ny=N;
92: // Initialize B-spline space
93: DMCreate(PETSC_COMM_WORLD,&user.iga);
94: DMSetType(user.iga, DMIGA);
95: DMIGAInitializeUniform2d(user.iga,PETSC_FALSE,2,3,
96: p,Nx,C,0.0,1.0,PETSC_TRUE,ng,
97: p,Ny,C,0.0,1.0,PETSC_TRUE,ng);
99: DMCreateGlobalVector(user.iga,&U);
100: FormInitialCondition(&user,U);
101: DMIGASetFieldName(user.iga, 0, "density");
102: DMIGASetFieldName(user.iga, 1, "velocity-u");
103: DMIGASetFieldName(user.iga, 2, "velocity-v");
105: DMGetMatrix(user.iga, MATAIJ, &J);
107: TSCreate(PETSC_COMM_WORLD,&ts);
108: TSSetType(ts,TSALPHA);
109: TSAlphaSetRadius(ts,0.5);
110: TSSetDM(ts,user.iga);
111: TSSetSolution(ts,U);
112: TSSetProblemType(ts,TS_NONLINEAR);
113: TSSetIFunction(ts,PETSC_NULL,FormResidual,&user);
114: TSSetIJacobian(ts,J,J,FormTangent,&user);
115: TSSetDuration(ts,1000000,1000.0);
116: TSSetInitialTimeStep(ts,0.0,0.001);
117: TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,PETSC_NULL);
118: TSMonitorSet(ts,OutputMonitor,PETSC_NULL,PETSC_NULL);
119: TSSetFromOptions(ts);
121: TSSolve(ts,U,&ftime);
122: TSGetTimeStepNumber(ts,&steps);
124: // Cleanup
125: TSDestroy(&ts);
126: DMDestroy(&user.iga);
127: PetscFinalize();
129: return 0;
130: }
134: PetscErrorCode FormInitialCondition(AppCtx *user,Vec U)
135: {
136: DMDALocalInfo info;
137: Field **u;
138: PetscScalar x,y,d1,d2,d3;
139: PetscInt i,j;
143: DMIGAGetLocalInfo(user->iga,&info);
144: DMIGAVecGetArray(user->iga,U,&u);
146: for(i=info.xs;i<info.xs+info.xm;i++){
147: x = user->L0*( (PetscScalar)i/(PetscScalar)info.mx );
148: for(j=info.ys;j<info.ys+info.ym;j++){
149: y = user->L0*( (PetscScalar)j/(PetscScalar)info.my );
151: d1 = sqrt(SQ(x-user->C1x)+SQ(y-user->C1y));
152: d2 = sqrt(SQ(x-user->C2x)+SQ(y-user->C2y));
153: d3 = sqrt(SQ(x-user->C3x)+SQ(y-user->C3y));
155: u[j][i].rho = -0.15+0.25*( tanh(0.5*(d1-user->R1)/user->Ca) +
156: tanh(0.5*(d2-user->R2)/user->Ca) +
157: tanh(0.5*(d3-user->R3)/user->Ca) );
158: u[j][i].ux = 0.0;
159: u[j][i].uy = 0.0;
160: }
161: }
162: DMIGAVecRestoreArray(user->iga,U,&u);
163: return(0);
164: }
168: PetscErrorCode FormResidual(TS ts,PetscReal t,Vec U,Vec Udot,Vec R,void *user)
169: {
174: DMDALocalInfo info;
175: DM dm;
176: Vec localU,localUdot,localR; // local versions
177: Field **h,**hdot,**r;
179: /* get the da from the snes */
180: TSGetDM(ts, &dm);
182: /* handle the vec U */
183: DMGetLocalVector(dm,&localU);
184: DMGlobalToLocalBegin(dm,U,INSERT_VALUES,localU);
185: DMGlobalToLocalEnd(dm,U,INSERT_VALUES,localU);
187: /* handle the vec Udot */
188: DMGetLocalVector(dm,&localUdot);
189: DMGlobalToLocalBegin(dm,Udot,INSERT_VALUES,localUdot);
190: DMGlobalToLocalEnd(dm,Udot,INSERT_VALUES,localUdot);
192: /* handle the vec R */
193: DMGetLocalVector(dm,&localR);
194: VecZeroEntries(localR);
196: /* Get the arrays from the vectors */
197: DMIGAVecGetArray(dm,localU,&h);
198: DMIGAVecGetArray(dm,localUdot,&hdot);
199: DMIGAVecGetArray(dm,localR,&r);
201: /* Grab the local info and call the local residual routine */
202: DMIGAGetLocalInfo(dm,&info);
203: FormResidualLocal(&info,t,h,hdot,r,(AppCtx *) user);
205: /* Restore the arrays */
206: DMIGAVecRestoreArray(dm,localR,&r);
207: DMIGAVecRestoreArray(dm,localUdot,&hdot);
208: DMIGAVecRestoreArray(dm,localU,&h);
210: /* Add contributions back to global R */
211: VecZeroEntries(R);
212: DMLocalToGlobalBegin(dm,localR,ADD_VALUES,R); // this one adds the values
213: DMLocalToGlobalEnd(dm,localR,ADD_VALUES,R);
215: /* Restore the local vectors */
216: DMRestoreLocalVector(dm,&localU);
217: DMRestoreLocalVector(dm,&localUdot);
218: DMRestoreLocalVector(dm,&localR);
220: return(0);
221: }
225: PetscErrorCode FormResidualLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,Field **r,AppCtx *user)
226: {
227: DM iga = user->iga;
228: BD bdX, bdY;
229: PetscInt px, py, ngx, ngy;
233: DMIGAGetPolynomialOrder(iga, &px, &py, PETSC_NULL);
234: DMIGAGetNumQuadraturePoints(iga, &ngx, &ngy, PETSC_NULL);
235: DMIGAGetBasisData(iga, &bdX, &bdY, PETSC_NULL);
237: // begin and end elements for this processor
238: PetscInt bex = bdX->own_b, eex = bdX->own_e;
239: PetscInt bey = bdY->own_b, eey = bdY->own_e;
241: // Allocate space for the 3D basis to be formed
242: PetscInt Nl = (px+1)*(py+1); // number of local basis functions
243: int numD = 6; // [0] = N, [1] = dN/dx, [2] = dN/dy
244: double **basis2D;
245: ierr= PetscMalloc(numD*sizeof(double*), &basis2D);
246: int i;
247: for(i=0;i<numD;i++) {
248: PetscMalloc(Nl*sizeof(double), &basis2D[i]);
249: }
251: PetscInt ind; // temporary index variable
252: PetscInt ie,je; // iterators for elements
253: PetscInt boffsetX,boffsetY; // offsets for l2g mapping
254: PetscInt ig,jg; // iterators for gauss points
255: PetscScalar gx,gy; // gauss point locations
256: PetscScalar wgtx,wgty,wgt; // gauss point weights
257: PetscInt iba,jba; // iterators for local basis (a, matrix rows)
258: PetscScalar Nx,dNx,dNxx,Ny,dNy,dNyy; // 1D basis functions and derivatives
259: PetscInt Ax,Ay; // global matrix row/col index
260: PetscScalar Na,Na_x,Na_y,Na_xx,Na_yy,Na_xy; // 2D basis for row loop (a)
262: PetscScalar R_rho,R_ux,R_uy;
263: PetscScalar rho,rho_t,rho_x,rho_y,rho_xx,rho_yy,rho_xy;
264: PetscScalar ux,ux_t,ux_x,ux_y;
265: PetscScalar uy,uy_t,uy_x,uy_y;
266: PetscScalar tau_xx,tau_xy,tau_yx,tau_yy;
267: PetscScalar p;
269: PetscScalar Ca2 = user->Ca*user->Ca;
270: PetscScalar rRe = 1.0/user->Re;
272: for(ie=bex;ie<=eex;ie++) { // Loop over elements
273: for(je=bey;je<=eey;je++) {
275: // get basis offsets used in the local-->global mapping
276: BDGetBasisOffset(bdX,ie,&boffsetX);
277: BDGetBasisOffset(bdY,je,&boffsetY);
279: for(ig=0;ig<ngx;ig++) { // Loop over gauss points
280: for(jg=0;jg<ngy;jg++) {
282: // Get gauss point locations and weights
283: // NOTE: gauss point and weight already mapped to the parameter space
284: BDGetGaussPt(bdX,ie,ig,&gx);
285: BDGetGaussPt(bdY,je,jg,&gy);
286: BDGetGaussWt(bdX,ie,ig,&wgtx);
287: BDGetGaussWt(bdY,je,jg,&wgty);
289: wgt = wgtx*wgty;
291: for(jba=0;jba<(py+1);jba++) { // Assemble the 2D basis
292: for(iba=0;iba<(px+1);iba++) {
294: BDGetBasis(bdX,ie,ig,iba,0,&Nx);
295: BDGetBasis(bdX,ie,ig,iba,1,&dNx);
296: BDGetBasis(bdX,ie,ig,iba,2,&dNxx);
298: BDGetBasis(bdY,je,jg,jba,0,&Ny);
299: BDGetBasis(bdY,je,jg,jba,1,&dNy);
300: BDGetBasis(bdY,je,jg,jba,2,&dNyy);
302: // 2D basis is a tensor product
303: ind = jba*(px+1)+iba;
304: basis2D[0][ind] = Nx * Ny; // N
305: basis2D[1][ind] = dNx * Ny; // dN/dx
306: basis2D[2][ind] = Nx * dNy; // dN/dy
307: basis2D[3][ind] = dNxx * Ny; // d^2N/dx^2
308: basis2D[4][ind] = Nx * dNyy; // d^2N/dy^2
309: basis2D[5][ind] = dNx * dNy; // d^2N/dxdy
311: }
312: } // end 2D basis assembly
314: // Problem coefficient evaluation
315: InterpolateSolution(basis2D,h,hdot,px,py,boffsetX,boffsetY,
316: &rho,&rho_t,&rho_x,&rho_y,
317: &rho_xx,&rho_yy,&rho_xy,
318: &ux,&ux_t,&ux_x,&ux_y,
319: &uy,&uy_t,&uy_x,&uy_y);
321: // compute pressure (Eq. 34.3)
322: p = 8.0/27.0*user->theta*rho/(1.0-rho)-rho*rho;
324: // compute viscous stress tensor (Eq. 34.4)
325: tau_xx = 2.0*ux_x - 2.0/3.0*(ux_x+uy_y);
326: tau_xy = ux_y + uy_x ;
327: tau_yy = 2.0*uy_y - 2.0/3.0*(ux_x+uy_y);
328: tau_yx = tau_xy;
330: for(jba=0;jba<(py+1);jba++) { // loop over basis 1st time (a, matrix rows)
331: for(iba=0;iba<(px+1);iba++) {
333: Ax = boffsetX+iba; // local to global map
334: Ay = boffsetY+jba;
336: ind = jba*(px+1)+iba;
337: Na = basis2D[0][ind];
338: Na_x = basis2D[1][ind];
339: Na_y = basis2D[2][ind];
340: Na_xx = basis2D[3][ind];
341: Na_yy = basis2D[4][ind];
342: Na_xy = basis2D[5][ind];
344: // (Eq. 19, modified to be dimensionless)
345: R_rho = Na*rho_t;
346: R_rho += -rho*(Na_x*ux + Na_y*uy);
348: R_ux = Na*ux*rho_t;
349: R_ux += Na*rho*ux_t;
350: R_ux += -rho*(Na_x*ux*ux + Na_y*ux*uy);
351: R_ux += -Na_x*p;
352: R_ux += rRe*(Na_x*tau_xx + Na_y*tau_xy);
353: //R_ux += -Ca2*rho*rho_x*(Na_xx+Na_xy);
354: //R_ux += -Ca2*Na_x*(rho_x*rho_x+rho_y*rho_y);
355: //R_ux += -Ca2*(rho_xx*Na + rho_x*Na_x + rho_xy*Na + rho_y*Na_x)*rho_x;
356: // rewritten uses Victor's corrections
357: R_ux += -Ca2*(Na_xx*rho_x + Na_xy*rho_y);
358: R_ux += -Ca2*Na_x*(rho_x*rho_x+rho_y*rho_y);
359: R_ux += -Ca2*Na*(rho_xx*rho_x+rho_xy*rho_y);
360: R_ux += -Ca2*rho_x*(Na_x*rho_x+Na_y*rho_y);
362: R_uy = Na*uy*rho_t;
363: R_uy += Na*rho*uy_t;
364: R_uy += -rho*(Na_x*uy*ux + Na_y*uy*uy);
365: R_uy += -Na_y*p;
366: R_uy += rRe*(Na_x*tau_yx + Na_y*tau_yy);
368: //R_uy += -Ca2*rho*rho_y*(Na_xy+Na_yy);
369: //R_uy += -Ca2*Na_y*(rho_x*rho_x+rho_y*rho_y);
370: //R_uy += -Ca2*(rho_xy*Na + rho_x*Na_y + rho_yy*Na + rho_y*Na_y)*rho_y;
372: R_uy += -Ca2*(Na_xy*rho_x + Na_yy*rho_y);
373: R_uy += -Ca2*Na_y*(rho_x*rho_x+rho_y*rho_y);
374: R_uy += -Ca2*Na*(rho_xy*rho_x+rho_yy*rho_y);
375: R_uy += -Ca2*rho_y*(Na_x*rho_x+Na_y*rho_y);
377: r[Ay][Ax].rho += R_rho*wgt;
378: r[Ay][Ax].ux += R_ux*wgt;
379: r[Ay][Ax].uy += R_uy*wgt;
381: }
382: } // end basis a loop
384: }
385: } // end gauss point loop
387: }
388: } // end element loop
390: for(i=0;i<numD;i++) {
391: PetscFree(basis2D[i]);
392: }
393: PetscFree(basis2D);
396: return(0);
397: }
401: PetscErrorCode FormTangent(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag,void *ctx)
402: {
406: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "FormTangent not implemented, use -snes_mf");
408: DMDALocalInfo info;
409: DM da_dof;
410: Vec localU,localUdot; // local versions
411: Field **h,**hdot;
413: /* get the da from the snes */
414: TSGetDM(ts,(DM*)&da_dof);
416: /* handle the vec U */
417: DMGetLocalVector(da_dof,&localU);
418: DMGlobalToLocalBegin(da_dof,U,INSERT_VALUES,localU);
419: DMGlobalToLocalEnd(da_dof,U,INSERT_VALUES,localU);
421: /* handle the vec Udot */
422: DMGetLocalVector(da_dof,&localUdot);
423: DMGlobalToLocalBegin(da_dof,Udot,INSERT_VALUES,localUdot);
424: DMGlobalToLocalEnd(da_dof,Udot,INSERT_VALUES,localUdot);
426: /* Get the arrays from the vectors */
427: DMDAVecGetArray(da_dof,localU,&h);
428: DMDAVecGetArray(da_dof,localUdot,&hdot);
430: /* Grab the local info and call the local tangent routine */
431: DMDAGetLocalInfo(da_dof,&info);
432: MatZeroEntries(*B); // pre-zero the matrix
433: FormTangentLocal(&info,t,h,hdot,shift,B,(AppCtx *) ctx);
434: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
435: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
436: if (*A != *B) { // then we could be matrix free
437: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
438: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
439: }
440: *flag = SAME_NONZERO_PATTERN; /* the sparsity pattern does not change */
442: DMDAVecRestoreArray(da_dof,localUdot,&hdot);
443: DMDAVecRestoreArray(da_dof,localU,&h);
445: /* Restore the arrays and local vectors */
446: DMRestoreLocalVector(da_dof,&localU);
447: DMRestoreLocalVector(da_dof,&localUdot);
449: return(0);
450: }
454: PetscErrorCode FormTangentLocal(DMDALocalInfo *info,PetscReal t,Field **h,Field **hdot,PetscReal shift,Mat *A,AppCtx *user)
455: {
458: return(0);
459: }
463: PetscErrorCode InterpolateSolution(double **basis2D,Field **x,Field **xdot,PetscInt px,PetscInt py,PetscInt boffsetX,PetscInt boffsetY,
464: PetscScalar *rho,PetscScalar *rho_t,PetscScalar *rho_x,PetscScalar *rho_y,
465: PetscScalar *rho_xx,PetscScalar *rho_yy,PetscScalar *rho_xy,
466: PetscScalar *ux,PetscScalar *ux_t,PetscScalar *ux_x,PetscScalar *ux_y,
467: PetscScalar *uy,PetscScalar *uy_t,PetscScalar *uy_x,PetscScalar *uy_y)
468: {
470: (*rho) = 0.0; (*rho_x) = 0.0; (*rho_y) = 0.0; (*rho_t) = 0.0;
471: (*rho_xx) = 0.0; (*rho_yy) = 0.0; (*rho_xy) = 0.0;
472: (*ux) = 0.0; (*ux_x) = 0.0; (*ux_y) = 0.0; (*ux_t) = 0.0;
473: (*uy) = 0.0; (*uy_x) = 0.0; (*uy_y) = 0.0; (*uy_t) = 0.0;
475: int ipa,jpa,ind;
476: for(jpa=0;jpa<(py+1);jpa++) {
477: for(ipa=0;ipa<(px+1);ipa++) {
479: ind = jpa*(px+1)+ipa;
480: (*rho) += basis2D[0][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
481: (*ux) += basis2D[0][ind] * x[boffsetY+jpa][boffsetX+ipa].ux;
482: (*uy) += basis2D[0][ind] * x[boffsetY+jpa][boffsetX+ipa].uy;
484: (*rho_x) += basis2D[1][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
485: (*ux_x) += basis2D[1][ind] * x[boffsetY+jpa][boffsetX+ipa].ux;
486: (*uy_x) += basis2D[1][ind] * x[boffsetY+jpa][boffsetX+ipa].uy;
488: (*rho_y) += basis2D[2][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
489: (*ux_y) += basis2D[2][ind] * x[boffsetY+jpa][boffsetX+ipa].ux;
490: (*uy_y) += basis2D[2][ind] * x[boffsetY+jpa][boffsetX+ipa].uy;
492: (*rho_xx) += basis2D[3][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
493: (*rho_yy) += basis2D[4][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
494: (*rho_xy) += basis2D[5][ind] * x[boffsetY+jpa][boffsetX+ipa].rho;
496: (*rho_t) += basis2D[0][ind] * xdot[boffsetY+jpa][boffsetX+ipa].rho;
497: (*ux_t) += basis2D[0][ind] * xdot[boffsetY+jpa][boffsetX+ipa].ux;
498: (*uy_t) += basis2D[0][ind] * xdot[boffsetY+jpa][boffsetX+ipa].uy;
500: }
501: }
504: return(0);
505: }
510: PetscErrorCode WriteSolution(Vec U, const char pattern[],int number)
511: {
513: PetscErrorCode ierr;
514: MPI_Comm comm;
515: char filename[256];
516: PetscViewer viewer;
519: sprintf(filename,pattern,number);
520: PetscObjectGetComm((PetscObject)U,&comm);
521: PetscViewerBinaryOpen(comm,filename,FILE_MODE_WRITE,&viewer);
522: VecView(U,viewer);
523: PetscViewerFlush(viewer);
524: PetscViewerDestroy(&viewer);
525: return(0);
526: }
530: PetscErrorCode OutputMonitor(TS ts,PetscInt it_number,PetscReal c_time,Vec U,void *mctx)
531: {
535: WriteSolution(U,"./nsk%d.dat",it_number);
537: return(0);
538: }