Actual source code: ex64.c

  1: static char help[] = "1D coupled Allen-Cahn and Cahn-Hilliard equation for constant mobility. Only c_v and eta are considered.\n\
  2: Runtime options include:\n\
  3: -xmin <xmin>\n\
  4: -xmax <xmax>\n\
  5: -ymin <ymin>\n\
  6: -T <T>, where <T> is the end time for the time domain simulation\n\
  7: -dt <dt>,where <dt> is the step size for the numerical integration\n\
  8: -gamma <gamma>\n\
  9: -theta_c <theta_c>\n\n";

 11: /*
 12: ./ex63 -ksp_type fgmres -snes_vi_monitor   -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason   -snes_ls_monitor -VG 10000000 -draw_fields 1,3,4  -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type svd  -ksp_monitor_true_residual -ksp_rtol 1.e-9 -snes_ls basic -T .0020 -P_casc .0005 -snes_monitor_solution -da_refine 10
 13:  */

 15:  #include petscsnes.h
 16:  #include petscdmda.h

 18: typedef struct{
 19:   PetscReal   dt,T; /* Time step and end time */
 20:   DM          da1,da1_clone,da2;
 21:   Mat         M;    /* Jacobian matrix */
 22:   Mat         M_0;
 23:   Vec         q,wv,cv,eta,DPsiv,DPsieta,logcv,logcv2;
 24:   Vec         work1,work2;
 25:   PetscScalar Mv,L,kaeta,kav,Evf,A,B,cv0,Sv;
 26:   PetscReal   xmin,xmax;
 27:   PetscInt    nx;
 28:   PetscBool   graphics;
 29:   PetscBool   periodic;
 30:   PetscBool   lumpedmass;
 31: }AppCtx;

 33: PetscErrorCode GetParams(AppCtx*);
 34: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
 35: PetscErrorCode SetUpMatrices(AppCtx*);
 36: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 37: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 38: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
 39: PetscErrorCode Update_q(AppCtx*);
 40: PetscErrorCode Update_u(Vec,AppCtx*);
 41: PetscErrorCode DPsi(AppCtx*);
 42: PetscErrorCode Llog(Vec,Vec);
 43: PetscErrorCode CheckRedundancy(SNES,IS,IS*,DM);

 47: int main(int argc, char **argv)
 48: {
 50:   Vec            x,r;  /* solution and residual vectors */
 51:   SNES           snes; /* Nonlinear solver context */
 52:   AppCtx         user; /* Application context */
 53:   Vec            xl,xu; /* Upper and lower bounds on variables */
 54:   Mat            J;
 55:   PetscScalar    t=0.0;
 56:   PetscViewer    view_out, view_p, view_q, view_psi, view_mat;
 57:   PetscReal      bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0};


 60:   PetscInitialize(&argc,&argv, (char*)0, help);
 61: 
 62:   /* Get physics and time parameters */
 63:   GetParams(&user);

 65:   if (user.periodic) {
 66:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 3, 1,PETSC_NULL,&user.da1);
 67:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 3, 1,PETSC_NULL,&user.da1_clone);
 68:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 1, 1,PETSC_NULL,&user.da2);
 69:   } else {
 70:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, -4, 3, 1,PETSC_NULL,&user.da1);
 71:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, -4, 3, 1,PETSC_NULL,&user.da1_clone);
 72:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, -4, 1, 1,PETSC_NULL,&user.da2);

 74:   }
 75:   /* Set Element type (triangular) */
 76:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
 77:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
 78: 
 79:   /* Set x and y coordinates */
 80:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
 81:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
 82:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
 83:   DMCreateGlobalVector(user.da1,&x);
 84:   VecDuplicate(x,&r);
 85:   VecDuplicate(x,&xl);
 86:   VecDuplicate(x,&xu);
 87:   VecDuplicate(x,&user.q);
 88: 
 89:   /* Get global vector user->wv from da2 and duplicate other vectors */
 90:   DMCreateGlobalVector(user.da2,&user.wv);
 91:   VecDuplicate(user.wv,&user.cv);
 92:   VecDuplicate(user.wv,&user.eta);
 93:   VecDuplicate(user.wv,&user.DPsiv);
 94:   VecDuplicate(user.wv,&user.DPsieta);
 95:   VecDuplicate(user.wv,&user.logcv);
 96:   VecDuplicate(user.wv,&user.logcv2);
 97:   VecDuplicate(user.wv,&user.work1);
 98:   VecDuplicate(user.wv,&user.work2);

100:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
101:   DMGetMatrix(user.da1,MATAIJ,&user.M);
102:   /* Get the (usual) mass matrix structure from da2 */
103:   DMGetMatrix(user.da2,MATAIJ,&user.M_0);
104:   /* Form the jacobian matrix and M_0 */
105:   SetUpMatrices(&user);
106:   SetInitialGuess(x,&user);
107:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
108: 
109:   SNESCreate(PETSC_COMM_WORLD,&snes);
110:   SNESSetDM(snes,user.da1);

112:   SNESSetFunction(snes,r,FormFunction,(void*)&user);
113:   SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

115:   SNESSetType(snes,SNESVI);
116:   SNESSetFromOptions(snes);
117:   SetVariableBounds(user.da1,xl,xu);
118:   //SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone);
119:   SNESVISetVariableBounds(snes,xl,xu);
120:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
121:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
122:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
123:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
124:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
125:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */

127:   if (user.graphics) {
128:     PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);

130:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
131:   }
132:   while(t<user.T) {

134:     char         filename[PETSC_MAX_PATH_LEN];
135:     PetscScalar  a = 1.0;
136:     PetscInt     i;
137:     PetscViewer  view;


140:     SNESSetFunction(snes,r,FormFunction,(void*)&user);
141:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

143:     DPsi(&user);
144:     VecView(user.DPsiv,view_psi);
145:     VecView(user.DPsieta,view_psi);

147:     Update_q(&user);
148:     VecView(user.q,view_q);
149:     MatView(user.M,view_mat);

151:     SNESSolve(snes,PETSC_NULL,x);
152:     VecView(x,view_out);


155:     if (user.graphics) {
156:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
157:     }
158: 
159:     PetscInt its;
160:     SNESGetIterationNumber(snes,&its);
161:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);

163:     Update_u(x,&user);

165:     for (i=0; i < (int)(user.T/a) ; i++) {
166:       if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
167:         sprintf(filename,"output_%f",t);
168:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
169:         VecView(user.cv,view);
170:         VecView(user.eta,view);
171:         PetscViewerDestroy(&view);
172:       }
173: 
174:     }

176:     t = t + user.dt;
177: 
178:   }

180: 
181:   PetscViewerDestroy(&view_out);
182:   PetscViewerDestroy(&view_p);
183:   PetscViewerDestroy(&view_q);
184:   PetscViewerDestroy(&view_psi);
185:   PetscViewerDestroy(&view_mat);
186:   VecDestroy(&x);
187:   VecDestroy(&r);
188:   VecDestroy(&xl);
189:   VecDestroy(&xu);
190:   VecDestroy(&user.q);
191:   VecDestroy(&user.wv);
192:   VecDestroy(&user.cv);
193:   VecDestroy(&user.eta);
194:   VecDestroy(&user.DPsiv);
195:   VecDestroy(&user.DPsieta);
196:   VecDestroy(&user.logcv);
197:   VecDestroy(&user.logcv2);
198:   VecDestroy(&user.work1);
199:   VecDestroy(&user.work2);
200:   MatDestroy(&user.M);
201:   MatDestroy(&user.M_0);
202:   DMDestroy(&user.da1);
203:   DMDestroy(&user.da1_clone);
204:   DMDestroy(&user.da2);
205:   SNESDestroy(&snes);
206:   PetscFinalize();
207:   return 0;
208: }

212: PetscErrorCode Update_u(Vec X,AppCtx *user)
213: {
215:   PetscInt       i,n;
216:   PetscScalar    *xx,*wv_p,*cv_p,*eta_p;
217: 
219:   VecGetLocalSize(user->wv,&n);
220: 
221:   VecGetArray(X,&xx);
222:   VecGetArray(user->wv,&wv_p);
223:   VecGetArray(user->cv,&cv_p);
224:   VecGetArray(user->eta,&eta_p);
225: 
226: 
227:   for(i=0;i<n;i++) {
228:     wv_p[i] = xx[3*i];
229:     cv_p[i] = xx[3*i+1];
230:     eta_p[i] = xx[3*i+2];
231:   }
232:   VecRestoreArray(X,&xx);
233:   VecRestoreArray(user->wv,&wv_p);
234:   VecRestoreArray(user->cv,&cv_p);
235:   VecRestoreArray(user->eta,&eta_p);
236: 
237:   return(0);
238: }

242: PetscErrorCode Update_q(AppCtx *user)
243: {
245:   PetscScalar    *q_p, *w1, *w2;
246:   PetscInt       n,i;

249: 
250:   VecGetArray(user->q,&q_p);
251:   VecGetArray(user->work1,&w1);
252:   VecGetArray(user->work2,&w2);
253:   VecGetLocalSize(user->wv,&n);

255:   MatMult(user->M_0,user->cv,user->work1);
256:   VecScale(user->work1,-1.0);
257:   for (i=0;i<n;i++) {
258:     q_p[3*i]=w1[i];
259:   }

261:   MatMult(user->M_0,user->DPsiv,user->work1);
262:   for (i=0;i<n;i++) {
263:     q_p[3*i+1]=w1[i];
264:   }

266:   VecCopy(user->DPsieta,user->work1);
267:   VecScale(user->work1,user->L*user->dt);
268:   VecAXPY(user->work1,-1.0,user->eta);
269:   MatMult(user->M_0,user->work1,user->work2);
270:   for (i=0;i<n;i++) {
271:     q_p[3*i+2]=w2[i];
272:   }

274:   VecRestoreArray(user->q,&q_p);
275:   VecRestoreArray(user->work1,&w1);
276:   VecRestoreArray(user->work2,&w2);
277: 
278:   return(0);
279: }

283: PetscErrorCode DPsi(AppCtx* user)
284: {
285:   PetscErrorCode  ierr;
286:   PetscScalar     Evf=user->Evf,A=user->A,B=user->B,cv0=user->cv0;
287:   PetscScalar     *cv_p,*eta_p,*logcv_p,*logcv2_p,*DPsiv_p,*DPsieta_p;
288:   PetscInt        n,i;


292:   VecGetLocalSize(user->cv,&n);
293:   VecGetArray(user->cv,&cv_p);
294:   VecGetArray(user->eta,&eta_p);
295:   VecGetArray(user->logcv,&logcv_p);
296:   VecGetArray(user->logcv2,&logcv2_p);
297:   VecGetArray(user->DPsiv,&DPsiv_p);
298:   VecGetArray(user->DPsieta,&DPsieta_p);

300:   Llog(user->cv,user->logcv);

302:   VecCopy(user->cv,user->work1);
303:   VecScale(user->work1,-1.0);
304:   VecShift(user->work1,1.0);
305:   Llog(user->work1,user->logcv2);

307:   for (i=0;i<n;i++)
308:   {
309:     DPsiv_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*(eta_p[i]+1.0)*(eta_p[i]+1.0)*( Evf + logcv_p[i] - logcv2_p[i]) - 2.0*A*(cv_p[i] - cv0)*eta_p[i]*(eta_p[i]+2.0)*(eta_p[i]-1.0)*(eta_p[i]-1.0) + 2.0*B*(cv_p[i] - 1.0)*eta_p[i]*eta_p[i];
310: 
311:     DPsieta_p[i] = 4.0*eta_p[i]*(eta_p[i]-1.0)*(eta_p[i]+1.0)*(Evf*cv_p[i] + cv_p[i]*logcv_p[i] + (1.0-cv_p[i])*logcv2_p[i] ) - A*(cv_p[i] - cv0)*(cv_p[i] - cv0)*(4.0*eta_p[i]*eta_p[i]*eta_p[i] - 6.0*eta_p[i] + 2.0) + 2.0*B*(cv_p[i]-1.0)*(cv_p[i]-1.0)*eta_p[i];
312: 
313:   }

315:   VecRestoreArray(user->cv,&cv_p);
316:   VecRestoreArray(user->eta,&eta_p);
317:   VecGetArray(user->logcv,&logcv_p);
318:   VecGetArray(user->logcv2,&logcv2_p);
319:   VecRestoreArray(user->DPsiv,&DPsiv_p);
320:   VecRestoreArray(user->DPsieta,&DPsieta_p);


323:   return(0);

325: }


330: PetscErrorCode Llog(Vec X, Vec Y)
331: {
332:   PetscErrorCode    ierr;
333:   PetscScalar       *x,*y;
334:   PetscInt          n,i;

337: 
338:   VecGetArray(X,&x);
339:   VecGetArray(Y,&y);
340:   VecGetLocalSize(X,&n);
341:   for (i=0;i<n;i++) {
342:     if (x[i] < 1.0e-12) {
343:       y[i] = log(1.0e-12);
344:     }
345:     else {
346:       y[i] = log(x[i]);
347:     }
348:   }
349: 
350:   return(0);
351: }

355: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
356: {
357:   PetscErrorCode    ierr;


360:   PetscInt         n,i,Mda;
361:   PetscScalar           *xx,*cv_p,*wv_p,*eta_p;
362:   PetscViewer      view_out;
363:   /* needed for the void growth case */
364:   PetscScalar       xmid,cv_v=1.0,cv_m=user->Sv*user->cv0,eta_v=1.0,eta_m=0.0,h,lambda;
365:   PetscInt          nele,nen,idx[2];
366:   const PetscInt    *ele;
367:   PetscScalar       x[2];
368:   Vec               coords;
369:   const PetscScalar *_coords;
370:   PetscScalar       xwidth = user->xmax - user->xmin;


374:   VecGetLocalSize(X,&n);
375: 
376: 
377:   DMDAGetInfo(user->da2,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
378:   DMDAGetGhostedCoordinates(user->da2,&coords);
379:   VecGetArrayRead(coords,&_coords);

381:   if (user->periodic) {
382:     h = (user->xmax-user->xmin)/Mda;
383:   } else {
384:     h = (user->xmax-user->xmin)/(Mda-1.0);
385:   }

387:   xmid = (user->xmax + user->xmin)/2.0;
388:   lambda = 4.0*h;
389: 
390:   DMDAGetElements(user->da2,&nele,&nen,&ele);
391:   for (i=0;i < nele; i++) {
392:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
393: 
394:     x[0] = _coords[idx[0]];
395:     x[1] = _coords[idx[1]];
396: 
397: 
398:     PetscInt k;
399:     PetscScalar vals_DDcv[2],vals_cv[2],vals_eta[2],s,hhr,r;
400:     for (k=0; k < 2 ; k++) {
401:       s = PetscAbs(x[k] - xmid);
402:       if (s <= xwidth*(5.0/64.0)) {
403:         vals_cv[k] = cv_v;
404:         vals_eta[k] = eta_v;
405:         vals_DDcv[k] = 0.0;
406:       } else if (s> xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
407:         //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
408:         r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
409:         hhr = 0.25*(-r*r*r + 3*r + 2);
410:         vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
411:         vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
412:         vals_DDcv[k] = (cv_v - cv_m)*r*6.0/(lambda*lambda);
413:       } else {
414:         vals_cv[k] = cv_m;
415:         vals_eta[k] = eta_m;
416:         vals_DDcv[k] = 0.0;
417:       }
418:     }
419: 
420:     VecSetValuesLocal(user->cv,2,idx,vals_cv,INSERT_VALUES);
421:     VecSetValuesLocal(user->eta,2,idx,vals_eta,INSERT_VALUES);
422:     VecSetValuesLocal(user->work2,2,idx,vals_DDcv,INSERT_VALUES);
423: 
424:     VecAssemblyBegin(user->cv);
425:     VecAssemblyEnd(user->cv);
426:     VecAssemblyBegin(user->eta);
427:     VecAssemblyEnd(user->eta);
428:     VecAssemblyBegin(user->work2);
429:     VecAssemblyEnd(user->work2);
430: 
431:     DMDARestoreElements(user->da2,&nele,&nen,&ele);
432:     VecRestoreArrayRead(coords,&_coords);
433:   }

435:   DPsi(user);
436:   VecCopy(user->DPsiv,user->wv);
437:   VecAXPY(user->wv,-2.0*user->kav,user->work2);

439:   VecGetArray(X,&xx);
440:   VecGetArray(user->wv,&wv_p);
441:   VecGetArray(user->cv,&cv_p);
442:   VecGetArray(user->eta,&eta_p);

444:   for (i=0;i<n/3;i++)
445:   {
446:     xx[3*i]=wv_p[i];
447:     xx[3*i+1]=cv_p[i];
448:     xx[3*i+2]=eta_p[i];
449:   }

451:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
452:   VecView(user->wv,view_out);
453:   VecView(user->cv,view_out);
454:   VecView(user->eta,view_out);
455:   PetscViewerDestroy(&view_out);

457:   VecRestoreArray(X,&xx);
458:   VecRestoreArray(user->wv,&wv_p);
459:   VecRestoreArray(user->cv,&cv_p);
460:   VecRestoreArray(user->eta,&eta_p);
461:   return(0);
462: }

466: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
467: {
469:   AppCtx         *user=(AppCtx*)ctx;
470: 
472:   MatMultAdd(user->M,X,user->q,F);
473:   return(0);
474: }

478: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
479: {
481:   AppCtx         *user=(AppCtx*)ctx;
482: 
484:   *flg = SAME_NONZERO_PATTERN;
485:   MatCopy(user->M,*J,*flg);
486:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
487:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
488:   return(0);
489: }
492: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
493: {
495:   PetscScalar    **l,**u;
496:   PetscInt       xs,xm;
497:   PetscInt       i;
498: 
500:   DMDAVecGetArrayDOF(da,xl,&l);
501:   DMDAVecGetArrayDOF(da,xu,&u);
502: 
503:   DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
504: 
505: 
506:   for(i=xs; i < xs+xm;i++) {
507:     l[i][0] = -SNES_VI_INF;
508:     l[i][1] = 0.0;
509:     l[i][2] = 0.0;
510:     u[i][0] = SNES_VI_INF;
511:     u[i][1] = 1.0;
512:     u[i][2] = 1.0;
513:   }
514: 
515:   DMDAVecRestoreArrayDOF(da,xl,&l);
516:   DMDAVecRestoreArrayDOF(da,xu,&u);
517:   return(0);
518: }

522: PetscErrorCode GetParams(AppCtx* user)
523: {
525:   PetscBool      flg;
526: 
528: 
529:   /* Set default parameters */
530:   user->xmin  = 0.0; user->xmax = 128.0;
531:   user->Mv    = 1.0;
532:   user->L     = 1.0;
533:   user->kaeta = 1.0;
534:   user->kav   = 0.5;
535:   user->Evf   = 9.09;
536:   user->A     = 9.09;
537:   user->B     = 9.09;
538:   user->cv0   = 1.13e-4;
539:   user->Sv    = 500.0;
540:   user->dt    = 1.0e-5;
541:   user->T     = 1.0e-2;
542:   user->graphics = PETSC_TRUE;
543:   user->periodic = PETSC_FALSE;
544:   user->lumpedmass = PETSC_FALSE;
545: 
546:   PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
547:   PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
548:   PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
549:   PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
550:   PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
551:   PetscOptionsBool("-periodic","Use periodic boundary conditions\n","None",user->periodic,&user->periodic,&flg);
552:   PetscOptionsBool("-lumpedmass","Use lumped mass matrix\n","None",user->lumpedmass,&user->lumpedmass,&flg);
553: 
554:   return(0);
555:  }


560: PetscErrorCode SetUpMatrices(AppCtx* user)
561: {
562:   PetscErrorCode    ierr;
563:   PetscInt          nele,nen,i,n;
564:   const PetscInt    *ele;
565:   PetscScalar       dt=user->dt,h;
566: 
567:   PetscInt          idx[2];
568:   PetscScalar       eM_0[2][2],eM_2[2][2];
569:   Mat               M=user->M;
570:   Mat               M_0=user->M_0;
571:   PetscInt          Mda;

573: 


577:   MatGetLocalSize(M,&n,PETSC_NULL);
578:   DMDAGetInfo(user->da1,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);

580:   if (user->periodic) {
581:     h = (user->xmax-user->xmin)/Mda;
582:   } else {
583:     h = (user->xmax-user->xmin)/(Mda-1.0);
584:   }
585:   if (user->lumpedmass) {
586:     eM_0[0][0] = h/2.0;
587:     eM_0[1][1] = h/2.0;
588:     eM_0[0][1] = eM_0[1][0] = 0.0;
589:   } else {
590:     eM_0[0][0]=eM_0[1][1]=h/3.0;
591:     eM_0[0][1]=eM_0[1][0]=h/6.0;
592:   }
593:   eM_2[0][0]=eM_2[1][1]=1.0/h;
594:   eM_2[0][1]=eM_2[1][0]=-1.0/h;

596:   /* Get local element info */
597:   DMDAGetElements(user->da1,&nele,&nen,&ele);
598:   for(i=0;i < nele;i++) {
599: 
600:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];


603:     PetscInt    row,cols[4],r,row_M_0,cols2[2];
604:     PetscScalar vals[4],vals_M_0[2],vals2[2];

606:     for(r=0;r<2;r++) {
607:       row_M_0 = idx[r];
608:       vals_M_0[0]=eM_0[r][0];
609:       vals_M_0[1]=eM_0[r][1];
610:       MatSetValuesLocal(M_0,1,&row_M_0,2,idx,vals_M_0,ADD_VALUES);

612:       row = 3*idx[r];
613:       cols[0] = 3*idx[0];     vals[0] = dt*eM_2[r][0]*user->Mv;
614:       cols[1] = 3*idx[1];     vals[1] = dt*eM_2[r][1]*user->Mv;
615:       cols[2] = 3*idx[0]+1;   vals[2] = eM_0[r][0];
616:       cols[3] = 3*idx[1]+1;   vals[3] = eM_0[r][1];
617: 
618:       /* Insert values in matrix M for 1st dof */
619:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
620: 
621:       row = 3*idx[r]+1;
622:       cols[0] = 3*idx[0];     vals[0] = -eM_0[r][0];
623:       cols[1] = 3*idx[1];     vals[1] = -eM_0[r][1];
624:       cols[2] = 3*idx[0]+1;   vals[2] = 2.0*user->kav*eM_2[r][0];
625:       cols[3] = 3*idx[1]+1;   vals[3] = 2.0*user->kav*eM_2[r][1];

627:       /* Insert values in matrix M for 2nd dof */
628:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
629: 
630: 
631:       row = 3*idx[r]+2;
632:       cols2[0] = 3*idx[0]+2;   vals2[0] = eM_0[r][0] + user->dt*2.0*user->L*user->kaeta*eM_2[r][0];
633:       cols2[1] = 3*idx[1]+2;   vals2[1] = eM_0[r][1] + user->dt*2.0*user->L*user->kaeta*eM_2[r][1];
634:       MatSetValuesLocal(M,1,&row,2,cols2,vals2,ADD_VALUES);
635:     }
636:   }

638:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
639:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
640: 
641:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
642:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
643: 
644:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
645: 
646: 
647:   return(0);
648: }

652: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
653: {
655:   PetscScalar    **uin,**uout;
656:   Vec            UIN, UOUT;
657:   PetscInt       xs,xm,*outindex;
658:   const PetscInt *index;
659:   PetscInt       k,i,l,n,M,cnt=0;
660: 
662:   DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
663:   DMGetGlobalVector(da,&UIN);
664:   VecSet(UIN,0.0);
665:   DMGetLocalVector(da,&UOUT);
666:   DMDAVecGetArrayDOF(da,UIN,&uin);
667:   ISGetIndices(act,&index);
668:   ISGetLocalSize(act,&n);
669:   for (k=0;k<n;k++){
670:     l = index[k]%5;
671:     i = index[k]/5;
672:     uin[i][l]=1.0;
673:   }
674:   printf("Number of active constraints before applying redundancy %d\n",n);
675:   ISRestoreIndices(act,&index);
676:   DMDAVecRestoreArrayDOF(da,UIN,&uin);
677:   DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
678:   DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
679:   DMDAVecGetArrayDOF(da,UOUT,&uout);
680: 
681:   DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

683:   for(i=xs; i < xs+xm;i++) {
684:     if (uout[i-1][1] && uout[i][1] && uout[i+1][1])
685:              uout[i][0] = 1.0;
686:     if (uout[i-1][3] && uout[i][3] && uout[i+1][3])
687:              uout[i][2] = 1.0;
688:   }

690:   for(i=xs; i < xs+xm;i++) {
691:     for(l=0;l<5;l++) {
692:       if (uout[i][l])
693:         cnt++;
694:     }
695:   }

697:   printf("Number of active constraints after applying redundancy %d\n",cnt);
698: 

700:   PetscMalloc(cnt*sizeof(PetscInt),&outindex);
701:   cnt = 0;
702: 
703:   for(i=xs; i < xs+xm;i++) {
704:     for(l=0;l<5;l++) {
705:       if (uout[i][l])
706:         outindex[cnt++] = 5*(i)+l;
707:     }
708:   }
709: 
710: 
711:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
712:   DMDAVecRestoreArrayDOF(da,UOUT,&uout);
713:   DMRestoreGlobalVector(da,&UIN);
714:   DMRestoreLocalVector(da,&UOUT);
715:   return(0);
716: }