Actual source code: ex60.c

  1: static char help[] = "2D coupled Allen-Cahn and Cahn-Hilliard equation for constant mobility and triangular elements.\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:  ./ex60 -ksp_type fgmres -pc_type mg  -snes_vi_monitor   -snes_atol 1.e-11  -da_grid_x 72 -da_grid_y 72 -ksp_rtol 1.e-8  -T 0.1  -VG 100 -pc_type lu -ksp_monitor_true_residual -pc_factor_mat_solver_package superlu -snes_converged_reason -ksp_converged_reason -pc_type sor  -ksp_rtol 1.e-9  -snes_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_monitor_solution

 14:  */

 16: /*
 17:    Possible additions to the code. At each iteration count the number of solution elements that are at the upper bound and stop the program if large

 19:    Add command-line option for constant or degenerate mobility 
 20:    Add command-line option for graphics at each time step

 22:    Check time-step business; what should it be? How to check that it is good? 
 23:    Make random right hand side forcing function proportional to time step so smaller time steps don't mean more radiation
 24:    How does the multigrid linear solver work now?
 25:    What happens when running with degenerate mobility

 27:    
 28:  */

 30:  #include petscsnes.h
 31:  #include petscdmda.h

 33: typedef struct{
 34:   PetscReal   dt,T; /* Time step and end time */
 35:   DM          da1,da2;
 36:   Mat         M;    /* Jacobian matrix */
 37:   Mat         M_0;
 38:   Vec         q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Rr,Riv;
 39:   Vec         work1,work2,work3,work4;
 40:   PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,P_casc,VG; /* physics parameters */
 41:   PetscReal   xmin,xmax,ymin,ymax;
 42:   PetscInt    Mda, Nda;
 43: }AppCtx;

 45: PetscErrorCode GetParams(AppCtx*);
 46: PetscErrorCode SetRandomVectors(AppCtx*);
 47: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
 48: PetscErrorCode SetUpMatrices(AppCtx*);
 49: PetscErrorCode UpdateMatrices(AppCtx*);
 50: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 51: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 52: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
 53: PetscErrorCode Update_q(AppCtx*);
 54: PetscErrorCode Update_u(Vec,AppCtx*);
 55: PetscErrorCode DPsi(AppCtx*);
 56: PetscErrorCode LaplacianFiniteDifference(AppCtx*);
 57: PetscErrorCode Llog(Vec,Vec);
 60: int main(int argc, char **argv)
 61: {
 63:   Vec            x,r;  /* Solution and residual vectors */
 64:   SNES           snes; /* Nonlinear solver context */
 65:   AppCtx         user; /* Application context */
 66:   Vec            xl,xu; /* Upper and lower bounds on variables */
 67:   Mat            J;
 68:   PetscScalar    t=0.0;
 69:   PetscViewer    view_out, view_q, view_psi, view_mat;
 70:   PetscViewer    view_rand;
 71:   IS             inactiveconstraints;
 72:   PetscInt       ninactiveconstraints,N;

 74:   PetscInitialize(&argc,&argv, (char*)0, help);
 75: 
 76:   /* Get physics and time parameters */
 77:   GetParams(&user);
 78:   /* Create a 1D DA with dof = 5; the whole thing */
 79:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 5, 1,PETSC_NULL,PETSC_NULL,&user.da1);
 80: 
 81:   /* Create a 1D DA with dof = 1; for individual componentes */
 82:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 1, 1,PETSC_NULL,PETSC_NULL,&user.da2);


 85:   /* Set Element type (triangular) */
 86:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
 87:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
 88: 
 89:   /* Set x and y coordinates */
 90:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
 91:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
 92:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
 93:   DMCreateGlobalVector(user.da1,&x);
 94:   VecGetSize(x,&N);
 95:   VecDuplicate(x,&r);
 96:   VecDuplicate(x,&xl);
 97:   VecDuplicate(x,&xu);
 98:   VecDuplicate(x,&user.q);
 99: 
100:   /* Get global vector user->wv from da2 and duplicate other vectors */
101:   DMCreateGlobalVector(user.da2,&user.wv);
102:   VecDuplicate(user.wv,&user.cv);
103:   VecDuplicate(user.wv,&user.wi);
104:   VecDuplicate(user.wv,&user.ci);
105:   VecDuplicate(user.wv,&user.eta);
106:   VecDuplicate(user.wv,&user.cvi);
107:   VecDuplicate(user.wv,&user.DPsiv);
108:   VecDuplicate(user.wv,&user.DPsii);
109:   VecDuplicate(user.wv,&user.DPsieta);
110:   VecDuplicate(user.wv,&user.Pv);
111:   VecDuplicate(user.wv,&user.Pi);
112:   VecDuplicate(user.wv,&user.Piv);
113:   VecDuplicate(user.wv,&user.logcv);
114:   VecDuplicate(user.wv,&user.logci);
115:   VecDuplicate(user.wv,&user.logcvi);
116:   VecDuplicate(user.wv,&user.work1);
117:   VecDuplicate(user.wv,&user.work2);
118:   VecDuplicate(user.wv,&user.Rr);
119:   VecDuplicate(user.wv,&user.Riv);
120: 

122:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
123:   DMGetMatrix(user.da1,MATAIJ,&user.M);
124:   /* Get the (usual) mass matrix structure from da2 */
125:   DMGetMatrix(user.da2,MATAIJ,&user.M_0);
126:   SetInitialGuess(x,&user);
127:   /* Form the jacobian matrix and M_0 */
128:   SetUpMatrices(&user);
129:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
130: 
131:   SNESCreate(PETSC_COMM_WORLD,&snes);
132:   SNESSetDM(snes,user.da1);

134:   SNESSetFunction(snes,r,FormFunction,(void*)&user);
135:   SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
136: 

138:    SNESSetType(snes,SNESVI);

140:   SNESSetFromOptions(snes);
141:   SetVariableBounds(user.da1,xl,xu);
142:   SNESVISetVariableBounds(snes,xl,xu);
143: 
144:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_rand",FILE_MODE_WRITE,&view_rand);
145:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
146:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
147:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
148:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
149: 
150:   while (t<user.T) {
151:     SNESSetFunction(snes,r,FormFunction,(void*)&user);
152:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

154:     SetRandomVectors(&user);
155:     /*    VecView(user.Pv,view_rand);
156:     VecView(user.Pi,view_rand);
157:      VecView(user.Piv,view_rand);*/

159:     DPsi(&user);
160:     /*    VecView(user.DPsiv,view_psi);
161:     VecView(user.DPsii,view_psi);
162:      VecView(user.DPsieta,view_psi);*/

164:     Update_q(&user);
165:     /*    VecView(user.q,view_q);
166:      MatView(user.M,view_mat);*/
167:     SNESSolve(snes,PETSC_NULL,x);
168:     SNESVIGetInactiveSet(snes,&inactiveconstraints);
169:     ISGetSize(inactiveconstraints,&ninactiveconstraints);
170:     /* if (ninactiveconstraints < .90*N) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,"To many active constraints, model has become non-physical"); */

172:     /*    VecView(x,view_out);*/
173:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
174:     PetscInt its;
175:     SNESGetIterationNumber(snes,&its);
176:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4f in %d iterations\n",t,its);

178:     Update_u(x,&user);
179:     UpdateMatrices(&user);
180:     t = t + user.dt;
181:   }
182: 
183:   PetscViewerDestroy(&view_rand);
184:   PetscViewerDestroy(&view_mat);
185:   PetscViewerDestroy(&view_q);
186:   PetscViewerDestroy(&view_out);
187:   PetscViewerDestroy(&view_psi);

189:   VecDestroy(&x);
190:   VecDestroy(&r);
191:   VecDestroy(&xl);
192:   VecDestroy(&xu);
193:   VecDestroy(&user.q);
194:   VecDestroy(&user.wv);
195:   VecDestroy(&user.cv);
196:   VecDestroy(&user.wi);
197:   VecDestroy(&user.ci);
198:   VecDestroy(&user.eta);
199:   VecDestroy(&user.cvi);
200:   VecDestroy(&user.DPsiv);
201:   VecDestroy(&user.DPsii);
202:   VecDestroy(&user.DPsieta);
203:   VecDestroy(&user.Pv);
204:   VecDestroy(&user.Pi);
205:   VecDestroy(&user.Piv);
206:   VecDestroy(&user.logcv);
207:   VecDestroy(&user.logci);
208:   VecDestroy(&user.logcvi);
209:   VecDestroy(&user.work1);
210:   VecDestroy(&user.work2);
211:   VecDestroy(&user.Rr);
212:   VecDestroy(&user.Riv);
213:   MatDestroy(&user.M);
214:   MatDestroy(&user.M_0);
215:   DMDestroy(&user.da1);
216:   DMDestroy(&user.da2);
217: 
218:   PetscFinalize();
219:   return 0;
220: }

224: PetscErrorCode Update_u(Vec X,AppCtx *user)
225: {
227:   PetscInt       i,n;
228:   PetscScalar    *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
229: 
231:   VecGetLocalSize(user->wv,&n);
232:   VecGetArray(X,&xx);
233:   VecGetArray(user->wv,&wv_p);
234:   VecGetArray(user->cv,&cv_p);
235:   VecGetArray(user->wi,&wi_p);
236:   VecGetArray(user->ci,&ci_p);
237:   VecGetArray(user->eta,&eta_p);


240:   for(i=0;i<n;i++) {
241:     wv_p[i] = xx[5*i];
242:     cv_p[i] = xx[5*i+1];
243:     wi_p[i] = xx[5*i+2];
244:     ci_p[i] = xx[5*i+3];
245:     eta_p[i] = xx[5*i+4];
246:   }
247:   VecRestoreArray(X,&xx);
248:   VecRestoreArray(user->wv,&wv_p);
249:   VecRestoreArray(user->cv,&cv_p);
250:   VecRestoreArray(user->wi,&wi_p);
251:   VecRestoreArray(user->ci,&ci_p);
252:   VecRestoreArray(user->eta,&eta_p);
253: 
254:   return(0);
255: }

259: PetscErrorCode Update_q(AppCtx *user)
260: {
262:   PetscScalar    *q_p,*w1,*w2;
263:   PetscInt       i,n;

266: 
267:   VecPointwiseMult(user->Rr,user->eta,user->eta);
268:   VecScale(user->Rr,user->Rsurf);
269:   VecShift(user->Rr,user->Rbulk);
270:   VecPointwiseMult(user->Riv,user->cv,user->ci);
271:   VecPointwiseMult(user->Riv,user->Rr,user->Riv);

273:   VecGetArray(user->q,&q_p);
274:   VecGetArray(user->work1,&w1);
275:   VecGetArray(user->work2,&w2);

277:   VecCopy(user->cv,user->work1);
278:   VecAXPY(user->work1,1.0,user->Pv);
279:   VecScale(user->work1,-1.0);
280:   MatMult(user->M_0,user->work1,user->work2);
281:   VecGetLocalSize(user->work1,&n);
282: 
283:  for (i=0;i<n;i++) {
284:        q_p[5*i]=w2[i];
285:   }
286: 
287:   MatMult(user->M_0,user->DPsiv,user->work1);
288:   for (i=0;i<n;i++) {
289:        q_p[5*i+1]=w1[i];
290:   }

292:   VecCopy(user->ci,user->work1);
293:   VecAXPY(user->work1,1.0,user->Pi);
294:   VecScale(user->work1,-1.0);
295:   MatMult(user->M_0,user->work1,user->work2);
296:   for (i=0;i<n;i++) {
297:        q_p[5*i+2]=w2[i];
298:   }

300:   MatMult(user->M_0,user->DPsii,user->work1);
301:   for (i=0;i<n;i++) {
302:        q_p[5*i+3]=w1[i];
303:   }

305:   VecCopy(user->eta,user->work1);
306:   VecScale(user->work1,-1.0/user->dt);
307:   VecAXPY(user->work1,user->L,user->DPsieta);
308:   VecAXPY(user->work1,-1.0,user->Piv);
309:   MatMult(user->M_0,user->work1,user->work2);
310:   for (i=0;i<n;i++) {
311:        q_p[5*i+4]=w2[i];
312:   }
313: 
314:   VecRestoreArray(user->q,&q_p);
315:   VecRestoreArray(user->work1,&w1);
316:   VecRestoreArray(user->work2,&w2);
317: 
318:   return(0);
319: }

323: PetscErrorCode DPsi(AppCtx* user)
324: {
325:   PetscErrorCode  ierr;
326:   PetscScalar     Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
327:   PetscScalar     *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
328:   PetscInt        n,i;


332:   VecGetLocalSize(user->cv,&n);
333:   VecGetArray(user->cv,&cv_p);
334:   VecGetArray(user->ci,&ci_p);
335:   VecGetArray(user->eta,&eta_p);
336:   VecGetArray(user->logcv,&logcv_p);
337:   VecGetArray(user->logci,&logci_p);
338:   VecGetArray(user->logcvi,&logcvi_p);
339:   VecGetArray(user->DPsiv,&DPsiv_p);
340:   VecGetArray(user->DPsii,&DPsii_p);
341:   VecGetArray(user->DPsieta,&DPsieta_p);
342: 
343:   Llog(user->cv,user->logcv);
344: 
345:   Llog(user->ci,user->logci);
346: 

348:   VecCopy(user->cv,user->cvi);
349:   VecAXPY(user->cvi,1.0,user->ci);
350:   VecScale(user->cvi,-1.0);
351:   VecShift(user->cvi,1.0);
352:   Llog(user->cvi,user->logcvi);
353: 
354:   for (i=0;i<n;i++)
355:   {
356:     DPsiv_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*( Evf + kBT*(logcv_p[i] - logcvi_p[i]) ) + eta_p[i]*eta_p[i]*2*A*(cv_p[i]-1);

358:     DPsii_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*( Eif + kBT*(logci_p[i] - logcvi_p[i]) ) + eta_p[i]*eta_p[i]*2*A*ci_p[i] ;
359: 
360:     DPsieta_p[i] = 2.0*(eta_p[i]-1.0)*( Evf*cv_p[i] + Eif*ci_p[i] + kBT*( cv_p[i]* logcv_p[i] + ci_p[i]* logci_p[i] + (1-cv_p[i]-ci_p[i])*logcvi_p[i] ) ) + 2.0*eta_p[i]*A*( (cv_p[i]-1.0)*(cv_p[i]-1.0) + ci_p[i]*ci_p[i]);
361: 
362: 
363:   }

365:   VecRestoreArray(user->cv,&cv_p);
366:   VecRestoreArray(user->ci,&ci_p);
367:   VecRestoreArray(user->eta,&eta_p);
368:   VecRestoreArray(user->logcv,&logcv_p);
369:   VecRestoreArray(user->logci,&logci_p);
370:   VecRestoreArray(user->logcvi,&logcvi_p);
371:   VecRestoreArray(user->DPsiv,&DPsiv_p);
372:   VecRestoreArray(user->DPsii,&DPsii_p);
373:   VecRestoreArray(user->DPsieta,&DPsieta_p);


376:   return(0);

378: }


383: PetscErrorCode Llog(Vec X, Vec Y)
384: {
385:   PetscErrorCode    ierr;
386:   PetscScalar       *x,*y;
387:   PetscInt          n,i;

390: 
391:   VecGetArray(X,&x);
392:   VecGetArray(Y,&y);
393:   VecGetLocalSize(X,&n);
394:   for (i=0;i<n;i++) {
395:     if (x[i] < 1.0e-12) {
396:       y[i] = log(1.0e-12);
397:     }
398:     else {
399:       y[i] = log(x[i]);
400:     }
401:   }
402: 
403:   return(0);
404: }


409: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
410: {
411:   PetscErrorCode    ierr;
412:   PetscInt          n,i;
413:   PetscScalar           *xx,*cv_p,*ci_p,*wv_p,*wi_p;
414:   PetscViewer       view;
415:   PetscScalar       initv = .00069;


419:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view);
420:   VecGetLocalSize(X,&n);

422:   VecSet(user->cv,initv);
423:   VecSet(user->ci,initv);
424:   VecSet(user->eta,0.0);

426:   DPsi(user);
427:   VecCopy(user->DPsiv,user->wv);
428:   VecCopy(user->DPsii,user->wi);

430:   VecGetArray(X,&xx);
431:   VecGetArray(user->cv,&cv_p);
432:   VecGetArray(user->ci,&ci_p);
433:   VecGetArray(user->wv,&wv_p);
434:   VecGetArray(user->wi,&wi_p);
435:   for (i=0;i<n/5;i++)
436:   {
437:     xx[5*i]=wv_p[i];
438:     xx[5*i+1]=cv_p[i];
439:     xx[5*i+2]=wi_p[i];
440:     xx[5*i+3]=ci_p[i];
441:     xx[5*i+4]=0.0;
442:   }

444:   VecView(user->wv,view);
445:   VecView(user->cv,view);
446:   VecView(user->wi,view);
447:   VecView(user->ci,view);
448:   PetscViewerDestroy(&view);

450:   VecRestoreArray(X,&xx);
451:   VecRestoreArray(user->cv,&cv_p);
452:   VecRestoreArray(user->ci,&ci_p);
453:   VecRestoreArray(user->wv,&wv_p);
454:   VecRestoreArray(user->wi,&wi_p);
455: 
456:   return(0);
457: }

461: PetscErrorCode SetRandomVectors(AppCtx* user)
462: {
464:   PetscInt       i,n,count=0;
465:   PetscScalar    *w1,*w2,*Pv_p,*eta_p;
466: 

469: 
470:   VecSetRandom(user->work1,PETSC_NULL);
471:   VecSetRandom(user->work2,PETSC_NULL);
472:   VecGetArray(user->work1,&w1);
473:   VecGetArray(user->work2,&w2);
474:   VecGetArray(user->Pv,&Pv_p);
475:   VecGetArray(user->eta,&eta_p);
476:   VecGetLocalSize(user->work1,&n);
477:   for (i=0;i<n;i++) {
478: 
479:     if (eta_p[i]>=0.8 || w1[i]>user->P_casc){
480:       Pv_p[i]=0;
481: 
482:     }
483:     else
484:     {
485:       Pv_p[i]=w2[i]*user->VG;
486:       count=count+1;
487:     }

489:   }
490: 
491:   VecCopy(user->Pv,user->Pi);
492:   VecScale(user->Pi,0.9);
493:   VecPointwiseMult(user->Piv,user->Pi,user->Pv);
494:   VecRestoreArray(user->work1,&w1);
495:   VecRestoreArray(user->work2,&w2);
496:   VecRestoreArray(user->Pv,&Pv_p);
497:   VecRestoreArray(user->eta,&eta_p);

499:   return(0);
500: 
501: }

505: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
506: {
508:   AppCtx         *user=(AppCtx*)ctx;
509: 
511:   MatMultAdd(user->M,X,user->q,F);
512:   return(0);
513: }

517: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
518: {
520:   AppCtx         *user=(AppCtx*)ctx;
521: 
523:   *flg = SAME_NONZERO_PATTERN;
524:   MatCopy(user->M,*J,*flg);
525:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
526:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
527:   return(0);
528: }
531: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
532: {
534:   PetscScalar    ***l,***u;
535:   PetscInt       xs,xm,ys,ym;
536:   PetscInt       j,i;
537: 
539:   DMDAVecGetArrayDOF(da,xl,&l);
540:   DMDAVecGetArrayDOF(da,xu,&u);
541: 
542:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
543: 
544:   for (j=ys; j<ys+ym; j++) {
545:     for(i=xs; i < xs+xm;i++) {
546:       l[j][i][0] = -SNES_VI_INF;
547:       l[j][i][1] = 0.0;
548:       l[j][i][2] = -SNES_VI_INF;
549:       l[j][i][3] = 0.0;
550:       l[j][i][4] = 0.0;
551:       u[j][i][0] = SNES_VI_INF;
552:       u[j][i][1] = 1.0;
553:       u[j][i][2] = SNES_VI_INF;
554:       u[j][i][3] = 1.0;
555:       u[j][i][4] = 1.0;
556:     }
557:   }
558: 
559: 
560:   DMDAVecRestoreArrayDOF(da,xl,&l);
561:   DMDAVecRestoreArrayDOF(da,xu,&u);
562:   return(0);
563: }

567: PetscErrorCode GetParams(AppCtx* user)
568: {
570:   PetscBool      flg;
571: 
573: 
574:   /* Set default parameters */
575:   user->xmin = 0.0; user->xmax = 1.0;
576:   user->ymin = 0.0; user->ymax = 1.0;
577:   user->Dv = 1.0; user->Di=4.0;
578:   user->Evf = 0.8; user->Eif = 1.2;
579:   user->A = 1.0;
580:   user->kBT = 0.11;
581:   user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
582:   user->Rsurf = 10.0; user->Rbulk = 1.0;
583:   user->L = 10.0; user->P_casc = 0.05;
584:   user->T = 1.0e-2;    user->dt = 1.0e-4;
585:   user->VG = 100.0;

587:   PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
588:   PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
589:   PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
590:   PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
591:   PetscOptionsGetReal(PETSC_NULL,"-P_casc",&user->P_casc,&flg);
592:   PetscOptionsGetReal(PETSC_NULL,"-VG",&user->VG,&flg);
593: 

595:   return(0);
596:  }


601: PetscErrorCode SetUpMatrices(AppCtx* user)
602: {
603:   PetscErrorCode    ierr;
604:   PetscInt          nele,nen,i,j,n;
605:   const PetscInt    *ele;
606:   PetscScalar       dt=user->dt,hx,hy;
607: 
608:   PetscInt          idx[3],*nodes, *connect, k;
609:   PetscScalar       eM_0[3][3],eM_2_even[3][3],eM_2_odd[3][3];
610:   PetscScalar       cv_sum, ci_sum;
611:   Mat               M=user->M;
612:   Mat               M_0=user->M_0;
613:   PetscInt          Mda=user->Mda, Nda=user->Nda, ld, rd, ru, lu;
614:   PetscScalar       *cv_p,*ci_p;
615: 
617: 
618:   /*  MatSetOption(M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
619:    MatSetOption(M_0,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);*/

621:   /* Create the mass matrix M_0 */
622:   VecGetArray(user->cv,&cv_p);
623:   VecGetArray(user->ci,&ci_p);
624:   MatGetLocalSize(M,&n,PETSC_NULL);
625:   DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);

627:   PetscMalloc((Mda+1)*(Nda+1)*sizeof(PetscInt),&nodes);
628:   PetscMalloc(Mda*Nda*2*3*sizeof(PetscInt),&connect);
629:   hx = (user->xmax-user->xmin)/Mda;
630:   hy = (user->ymax-user->ymin)/Nda;
631:   for (j=0;j < Nda;j++) {
632:     for (i=0;i < Mda;i++) {
633:       nodes[j*(Mda+1)+i] = j*Mda+i;
634:     }
635:     nodes[j*(Mda+1)+Mda] = j*Mda;
636:   }
637: 
638:   for (i=0;i < Mda;i++){
639:     nodes[Nda*(Mda+1)+i] = i;
640:   }

642:   nodes[Nda*(Mda+1)+Mda] = 0;
643: 
644: 
645:   k = 0;
646:   for (j=0;j<Nda;j++) {
647:     for (i=0;i<Mda;i++) {
648: 
649:       /* ld = nodes[j][i]; */
650:       ld = nodes[j*(Mda+1)+i];
651:       /* rd = nodes[j+1][i]; */
652:       rd = nodes[(j+1)*(Mda+1)+i];
653:       /* ru = nodes[j+1][i+1]; */
654:       ru = nodes[(j+1)*(Mda+1)+i+1];
655:       /* lu = nodes[j][i+1]; */
656:       lu = nodes[j*(Mda+1)+i+1];

658:       /* connect[k][0]=ld; */
659:       connect[k*6]=ld;
660:       /* connect[k][1]=lu; */
661:       connect[k*6+1]=lu;
662:       /* connect[k][2]=ru; */
663:       connect[k*6+2]=rd;
664:       connect[k*6+3]=lu;
665:       connect[k*6+4]=ru;
666:       connect[k*6+5]=rd;
667: 
668:       k = k+1;
669:     }
670:   }
671: 
672: 
673:   eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hy/12.0;
674:   eM_0[0][1]=eM_0[0][2]=eM_0[1][0]=eM_0[1][2]=eM_0[2][0]=eM_0[2][1]=hx*hy/24.0;
675: 
676:   eM_2_odd[0][0] = 1.0;
677:   eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
678:   eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
679:   eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;

681:   eM_2_even[1][1] = 1.0;
682:   eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
683:   eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
684:   eM_2_even[0][2] = eM_2_even[2][0] = 0.0;


687:   for(k=0;k < Mda*Nda*2;k++) {
688:     idx[0] = connect[k*3];
689:     idx[1] = connect[k*3+1];
690:     idx[2] = connect[k*3+2];
691: 
692:     PetscInt    row,cols[6],r,row_M_0,cols3[3];
693:     PetscScalar vals[6],vals_M_0[3],vals3[3];
694: 
695:     for(r=0;r<3;r++) {
696:       row_M_0 = connect[k*3+r];
697: 
698:       vals_M_0[0]=eM_0[r][0];
699:       vals_M_0[1]=eM_0[r][1];
700:       vals_M_0[2]=eM_0[r][2];
701: 
702: 
703:       MatSetValues(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
704: 
705:       //cv_sum = (cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
706:       //ci_sum = (ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
707:       cv_sum = .0000069*user->Dv/user->kBT;
708:       ci_sum = .0000069*user->Di/user->kBT;
709: 
710:       if (k%2 == 0)  {
711: 
712:         row = 5*idx[r];
713:         cols[0] = 5*idx[0];     vals[0] = dt*eM_2_odd[r][0]*cv_sum;
714:         cols[1] = 5*idx[1];     vals[1] = dt*eM_2_odd[r][1]*cv_sum;
715:         cols[2] = 5*idx[2];     vals[2] = dt*eM_2_odd[r][2]*cv_sum;
716:         cols[3] = 5*idx[0]+1;   vals[3] = eM_0[r][0];
717:         cols[4] = 5*idx[1]+1;   vals[4] = eM_0[r][1];
718:         cols[5] = 5*idx[2]+1;   vals[5] = eM_0[r][2];

720: 
721:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
722: 
723: 
724:         row = 5*idx[r]+1;
725:         cols[0] = 5*idx[0];     vals[0] = -1.0*eM_0[r][0];
726:         cols[1] = 5*idx[1];     vals[1] = -1.0*eM_0[r][1];
727:         cols[2] = 5*idx[2];     vals[2] = -1.0*eM_0[r][2];
728:         cols[3] = 5*idx[0]+1;   vals[3] =  user->kav*eM_2_odd[r][0];
729:         cols[4] = 5*idx[1]+1;   vals[4] =  user->kav*eM_2_odd[r][1];
730:         cols[5] = 5*idx[2]+1;   vals[5] =  user->kav*eM_2_odd[r][2];
731: 
732:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
733: 
734:         row = 5*idx[r]+2;
735:         cols[0] = 5*idx[0]+2;   vals[0] =  dt*eM_2_odd[r][0]*ci_sum;
736:         cols[1] = 5*idx[1]+2;   vals[1] =  dt*eM_2_odd[r][1]*ci_sum;
737:         cols[2] = 5*idx[2]+2;   vals[2] =  dt*eM_2_odd[r][2]*ci_sum;
738:         cols[3] = 5*idx[0]+3;   vals[3] =  eM_0[r][0];
739:         cols[4] = 5*idx[1]+3;   vals[4] =  eM_0[r][1];
740:         cols[5] = 5*idx[2]+3;   vals[5] =  eM_0[r][2];
741: 
742:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
743: 
744: 
745:         row = 5*idx[r]+3;
746:         cols[0] = 5*idx[0]+2;   vals[0] = -1.0*eM_0[r][0];
747:         cols[1] = 5*idx[1]+2;   vals[1] = -1.0*eM_0[r][1];
748:         cols[2] = 5*idx[2]+2;   vals[2] = -1.0*eM_0[r][2];
749:         cols[3] = 5*idx[0]+3;   vals[3] =  user->kai*eM_2_odd[r][0];
750:         cols[4] = 5*idx[1]+3;   vals[4] =  user->kai*eM_2_odd[r][1];
751:         cols[5] = 5*idx[2]+3;   vals[5] =  user->kai*eM_2_odd[r][2];
752: 
753:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
754: 
755: 
756:         row = 5*idx[r]+4;
757:         /*
758:         cols3[0] = 5*idx[0]+4;   vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_odd[r][0];
759:         cols3[1] = 5*idx[1]+4;   vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_odd[r][1];
760:         cols3[2] = 5*idx[2]+4;   vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_odd[r][2];
761:          */
762:         cols3[0] = 5*idx[0]+4;   vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2_odd[r][0];
763:         cols3[1] = 5*idx[1]+4;   vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2_odd[r][1];
764:         cols3[2] = 5*idx[2]+4;   vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*eM_2_odd[r][2];
765: 
766:         MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
767: 
768: 
769:       }
770: 
771:       else {
772: 
773: 
774:         row = 5*idx[r];
775:         cols[0] = 5*idx[0];     vals[0] = dt*eM_2_even[r][0]*cv_sum;
776:         cols[1] = 5*idx[1];     vals[1] = dt*eM_2_even[r][1]*cv_sum;
777:         cols[2] = 5*idx[2];     vals[2] = dt*eM_2_even[r][2]*cv_sum;
778:         cols[3] = 5*idx[0]+1;   vals[3] = eM_0[r][0];
779:         cols[4] = 5*idx[1]+1;   vals[4] = eM_0[r][1];
780:         cols[5] = 5*idx[2]+1;   vals[5] = eM_0[r][2];


783: 
784:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
785: 
786: 
787:         row = 5*idx[r]+1;
788:         cols[0] = 5*idx[0];     vals[0] = -1.0*eM_0[r][0];
789:         cols[1] = 5*idx[1];     vals[1] = -1.0*eM_0[r][1];
790:         cols[2] = 5*idx[2];     vals[2] = -1.0*eM_0[r][2];
791:         cols[3] = 5*idx[0]+1;   vals[3] =  user->kav*eM_2_even[r][0];
792:         cols[4] = 5*idx[1]+1;   vals[4] =  user->kav*eM_2_even[r][1];
793:         cols[5] = 5*idx[2]+1;   vals[5] =  user->kav*eM_2_even[r][2];
794: 
795:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
796: 
797:         row = 5*idx[r]+2;
798:         cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2_even[r][0]*ci_sum;
799:         cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2_even[r][1]*ci_sum;
800:         cols[2] = 5*idx[2]+2;   vals[2] = dt*eM_2_even[r][2]*ci_sum;
801:         cols[3] = 5*idx[0]+3;   vals[3] = eM_0[r][0];
802:         cols[4] = 5*idx[1]+3;   vals[4] = eM_0[r][1];
803:         cols[5] = 5*idx[2]+3;   vals[5] = eM_0[r][2];
804: 
805:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
806: 
807:         row = 5*idx[r]+3;
808:         cols[0] = 5*idx[0]+2;   vals[0] = -1.0*eM_0[r][0];
809:         cols[1] = 5*idx[1]+2;   vals[1] = -1.0*eM_0[r][1];
810:         cols[2] = 5*idx[2]+2;   vals[2] = -1.0*eM_0[r][2];
811:         cols[3] = 5*idx[0]+3;   vals[3] =  user->kai*eM_2_even[r][0];
812:         cols[4] = 5*idx[1]+3;   vals[4] =  user->kai*eM_2_even[r][1];
813:         cols[5] = 5*idx[2]+3;   vals[5] =  user->kai*eM_2_even[r][2];
814: 
815:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
816: 
817:         row = 5*idx[r]+4;
818:         /*
819:         cols3[0] = 5*idx[0]+4;   vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_even[r][0];
820:         cols3[1] = 5*idx[1]+4;   vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_even[r][1];
821:         cols3[2] = 5*idx[2]+4;   vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_even[r][2];
822:          */
823:         cols3[0] = 5*idx[0]+4;   vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2_even[r][0];
824:         cols3[1] = 5*idx[1]+4;   vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2_even[r][1];
825:         cols3[2] = 5*idx[2]+4;   vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*eM_2_even[r][2];
826:         MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
827: 
828:       }
829: 
830:     }
831:   }

833:   PetscFree(nodes);
834:   PetscFree(connect);

836:   VecRestoreArray(user->cv,&cv_p);
837:   VecRestoreArray(user->ci,&ci_p);
838:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
839:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
840: 
841:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
842:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
843: 
844:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
845: 

847:   return(0);
848: }


853: PetscErrorCode UpdateMatrices(AppCtx* user)
854: {
855:   PetscErrorCode    ierr;
856:   PetscInt          i,j,n,Mda,Nda;
857: 
858:   PetscInt          idx[3],*nodes,*connect,k;
859:   PetscInt          ld,rd,lu,ru;
860:   PetscScalar       eM_2_odd[3][3],eM_2_even[3][3],h,dt=user->dt;
861:   Mat               M=user->M;
862:   PetscScalar       *cv_p,*ci_p,cv_sum,ci_sum;
864: 
865:   /* Create the mass matrix M_0 */
866:   MatGetLocalSize(M,&n,PETSC_NULL);
867:   VecGetArray(user->cv,&cv_p);
868:   VecGetArray(user->ci,&ci_p);
869:   DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);

871:   PetscMalloc((Mda+1)*(Nda+1)*sizeof(PetscInt),&nodes);
872:   PetscMalloc(Mda*Nda*2*3*sizeof(PetscInt),&connect);
873: 
874:   h = 100.0/Mda;

876:   for (j=0;j < Nda;j++) {
877:     for (i=0;i < Mda;i++) {
878:       nodes[j*(Mda+1)+i] = j*Mda+i;
879:     }
880:     nodes[j*(Mda+1)+Mda] = j*Mda;
881:   }
882:   for (i=0;i < Mda;i++){
883:     nodes[Nda*(Mda+1)+i]=i;
884:   }
885:   nodes[Nda*(Mda+1)+Mda]=0;

887: 
888:   k = 0;
889:   for (j=0;j<Nda;j++) {
890:     for (i=0;i<Mda;i++) {
891:       ld = nodes[j*(Mda+1)+i];
892:       rd = nodes[(j+1)*(Mda+1)+i];
893:       ru = nodes[(j+1)*(Mda+1)+i+1];
894:       lu = nodes[j*(Mda+1)+i+1];
895:       connect[k*6]=ld;
896:       connect[k*6+1]=lu;
897:       connect[k*6+2]=rd;
898:       connect[k*6+3]=lu;
899:       connect[k*6+4]=ru;
900:       connect[k*6+5]=rd;
901:       k = k+1;
902:     }
903:   }
904: 
905:   for(k=0;k < Mda*Nda*2;k++) {
906:     idx[0] = connect[k*3];
907:     idx[1] = connect[k*3+1];
908:     idx[2] = connect[k*3+2];
909: 
910:     PetscInt r,row,cols[3];
911:     PetscScalar vals[3];
912:     for (r=0;r<3;r++) {
913:       row = 5*idx[r];
914:       cols[0] = 5*idx[0];     vals[0] = 0.0;
915:       cols[1] = 5*idx[1];     vals[1] = 0.0;
916:       cols[2] = 5*idx[2];     vals[2] = 0.0;

918:       /* Insert values in matrix M for 1st dof */
919:       MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
920: 
921:       row = 5*idx[r]+2;
922:       cols[0] = 5*idx[0]+2;   vals[0] = 0.0;
923:       cols[1] = 5*idx[1]+2;   vals[1] = 0.0;
924:       cols[2] = 5*idx[2]+2;   vals[2] = 0.0;

926:       /* Insert values in matrix M for 3nd dof */
927:       MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
928:     }
929:   }
930: 
931:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
932:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

934:   eM_2_odd[0][0] = 1.0;
935:   eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
936:   eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
937:   eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;

939:   eM_2_even[1][1] = 1.0;
940:   eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
941:   eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
942:   eM_2_even[0][2] = eM_2_even[2][0] = 0.0;

944: 
945:   /* Get local element info */
946:   for(k=0;k < Mda*Nda*2;k++) {
947:     idx[0] = connect[k*3];
948:       idx[1] = connect[k*3+1];
949:       idx[2] = connect[k*3+2];
950: 
951:       PetscInt    row,cols[3],r;
952:       PetscScalar vals[3];
953: 
954:       for(r=0;r<3;r++) {
955: 
956:         // cv_sum = (1.0e-3+cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
957:         //ci_sum = (1.0e-3+ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
958:         cv_sum = .0000069*user->Dv/(user->kBT);
959:         ci_sum = .0000069*user->Di/user->kBT;

961:             if (k%2 == 0) /* odd triangle */ {
962: 
963:                 row = 5*idx[r];
964:                 cols[0] = 5*idx[0];     vals[0] = dt*eM_2_odd[r][0]*cv_sum;
965:                 cols[1] = 5*idx[1];     vals[1] = dt*eM_2_odd[r][1]*cv_sum;
966:                 cols[2] = 5*idx[2];     vals[2] = dt*eM_2_odd[r][2]*cv_sum;
967: 
968:                 /* Insert values in matrix M for 1st dof */
969:                 MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
970: 
971: 
972:                 row = 5*idx[r]+2;
973:                 cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2_odd[r][0]*ci_sum;
974:                 cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2_odd[r][1]*ci_sum;
975:                 cols[2] = 5*idx[2]+2;   vals[2] = dt*eM_2_odd[r][2]*ci_sum;
976: 
977:                 MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
978: 
979:             }
980: 
981:             else {
982:                 row = 5*idx[r];
983:                 cols[0] = 5*idx[0];     vals[0] = dt*eM_2_even[r][0]*cv_sum;
984:                 cols[1] = 5*idx[1];     vals[1] = dt*eM_2_even[r][1]*cv_sum;
985:                 cols[2] = 5*idx[2];     vals[2] = dt*eM_2_even[r][2]*cv_sum;
986: 
987:                 /* Insert values in matrix M for 1st dof */
988:                 MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
989: 
990: 
991:                 row = 5*idx[r]+2;
992:                 cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2_even[r][0]*ci_sum;
993:                 cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2_even[r][1]*ci_sum;
994:                 cols[2] = 5*idx[2]+2;   vals[2] = dt*eM_2_even[r][2]*ci_sum;
995:                 /* Insert values in matrix M for 3nd dof */
996:                 MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
997: 
998:             }
999:         }
1000: 
1001:     }

1003:   PetscFree(nodes);
1004:   PetscFree(connect);
1005:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1006:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1007:   VecRestoreArray(user->cv,&cv_p);
1008:   VecRestoreArray(user->ci,&ci_p);


1011:   return(0);
1012: }