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: }