Actual source code: ex61.c

petsc-3.4.4 2014-03-13
  1: static char help[] = "2D coupled Allen-Cahn and Cahn-Hilliard equation for constant mobility and triangular elements. Use periodic boundary condidtions.\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:  ./ex61 -ksp_type gmres -snes_vi_monitor   -snes_atol 1.e-11  -da_refine 3  -T 0.1   -ksp_monitor_true_residual -pc_type lu -pc_factor_mat_solver_package superlu -snes_converged_reason -ksp_converged_reason  -ksp_rtol 1.e-9  -snes_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_linesearch_type basic

 14: ./ex61 -ksp_type gmres -snes_vi_monitor   -snes_atol 1.e-11  -da_refine 4 -T 0.1   -ksp_monitor_true_residual -pc_type sor -snes_converged_reason -ksp_converged_reason  -ksp_rtol 1.e-9  -snes_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_linesearch_type basic

 16: ./ex61 -ksp_type fgmres -snes_vi_monitor   -snes_atol 1.e-11  -da_refine 5 -T 0.1   -ksp_monitor_true_residual -snes_converged_reason -ksp_converged_reason  -ksp_rtol 1.e-9  -snes_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_linesearch_type basic -pc_type mg -pc_mg_galerkin

 18: ./ex61 -ksp_type fgmres -snes_vi_monitor   -snes_atol 1.e-11  -da_refine 5 -snes_converged_reason -ksp_converged_reason   -snes_linesearch_monitor -VG 1 -draw_fields 1,3,4  -pc_type mg -pc_mg_galerkin -log_summary -dt .0000000000001 -mg_coarse_pc_type svd  -ksp_monitor_true_residual -ksp_rtol 1.e-9

 20: Movie version
 21: ./ex61 -ksp_type fgmres -snes_vi_monitor   -snes_atol 1.e-11  -da_refine 6 -snes_converged_reason -ksp_converged_reason   -snes_linesearch_monitor -VG 10 -draw_fields 1,3,4  -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type redundant -mg_coarse_redundant_pc_type svd  -ksp_monitor_true_residual -ksp_rtol 1.e-9 -snes_linesearch_type basic -T .0020

 23:  */

 25: /*
 26:    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

 28:    Is the solution at time 0 nonsense?  Looks totally different from first time step. Why does cubic line search at beginning screw it up?

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

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


 39:  */

 41:  #include petscsnes.h
 42:  #include petscdmda.h

 44: typedef struct {
 45:   PetscReal   dt,T; /* Time step and end time */
 46:   PetscReal   dtevent;  /* time scale of radiation events, roughly one event per dtevent */
 47:   PetscInt    maxevents; /* once this number of events is reached no more events are generated */
 48:   PetscReal   initv;    /* initial value of phase variables */
 49:   PetscReal   initeta;
 50:   PetscBool   degenerate;  /* use degenerate mobility */
 51:   PetscReal   smallnumber;
 52:   PetscBool   graphics;
 53:   PetscInt    domain;
 54:   PetscBool   radiation;
 55:   PetscBool   voidgrowth; /* use initial conditions for void growth */
 56:   DM          da1,da2;
 57:   Mat         M;    /* Jacobian matrix */
 58:   Mat         M_0;
 59:   Vec         q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Riv;
 60:   Vec         phi1,phi2,Phi2D_V,Sv,Si; /* for twodomain modeling */
 61:   Vec         work1,work2,work3,work4;
 62:   PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,VG; /* physics parameters */
 63:   PetscScalar Svr,Sir,cv_eq,ci_eq; /* for twodomain modeling */
 64:   PetscReal   asmallnumber; /* gets added to degenerate mobility */
 65:   PetscReal   xmin,xmax,ymin,ymax;
 66:   PetscInt    Mda, Nda;
 67:   PetscViewer graphicsfile;  /* output of solution at each times step */
 68: } AppCtx;

 70: PetscErrorCode GetParams(AppCtx*);
 71: PetscErrorCode SetRandomVectors(AppCtx*,PetscReal);
 72: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
 73: PetscErrorCode SetUpMatrices(AppCtx*);
 74: PetscErrorCode UpdateMatrices(AppCtx*);
 75: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 76: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 77: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
 78: PetscErrorCode Update_q(AppCtx*);
 79: PetscErrorCode Update_u(Vec,AppCtx*);
 80: PetscErrorCode DPsi(AppCtx*);
 81: PetscErrorCode LaplacianFiniteDifference(AppCtx*);
 82: PetscErrorCode Llog(Vec,Vec);
 83: PetscErrorCode Phi(AppCtx*);
 84: PetscErrorCode Phi_read(AppCtx*);

 88: int main(int argc, char **argv)
 89: {
 90:   PetscErrorCode      ierr;
 91:   Vec                 x,r;  /* Solution and residual vectors */
 92:   SNES                snes; /* Nonlinear solver context */
 93:   AppCtx              user; /* Application context */
 94:   Vec                 xl,xu; /* Upper and lower bounds on variables */
 95:   Mat                 J;
 96:   PetscScalar         t=0.0,normq;
 97:   IS                  inactiveconstraints;
 98:   PetscInt            ninactiveconstraints,N;
 99:   SNESConvergedReason reason;
100:   char                cv_filename[80],eta_filename[80];
101:   /*PetscViewer         view_out, view_cv,view_eta,view_vtk_cv,view_vtk_eta;*/
102:   /*PetscReal           bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0}; */
103:   /*  PetscViewer         view_out, view_q, view_psi, view_mat;*/
104:   /*  PetscViewer         view_rand;*/

106:   PetscInitialize(&argc,&argv, (char*)0, help);

108:   /* Get physics and time parameters */
109:   GetParams(&user);
110:   /* Create a 1D DA with dof = 5; the whole thing */
111:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,&user.da1);

113:   /* Create a 1D DA with dof = 1; for individual componentes */
114:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&user.da2);


117:   /* Set Element type (triangular) */
118:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
119:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);

121:   /* Set x and y coordinates */
122:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);
123:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);


126:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
127:   DMCreateGlobalVector(user.da1,&x);
128:   VecGetSize(x,&N);
129:   VecDuplicate(x,&r);
130:   VecDuplicate(x,&xl);
131:   VecDuplicate(x,&xu);
132:   VecDuplicate(x,&user.q);

134:   /* Get global vector user->wv from da2 and duplicate other vectors */
135:   DMCreateGlobalVector(user.da2,&user.wv);
136:   VecDuplicate(user.wv,&user.cv);
137:   VecDuplicate(user.wv,&user.wi);
138:   VecDuplicate(user.wv,&user.ci);
139:   VecDuplicate(user.wv,&user.eta);
140:   VecDuplicate(user.wv,&user.cvi);
141:   VecDuplicate(user.wv,&user.DPsiv);
142:   VecDuplicate(user.wv,&user.DPsii);
143:   VecDuplicate(user.wv,&user.DPsieta);
144:   VecDuplicate(user.wv,&user.Pv);
145:   VecDuplicate(user.wv,&user.Pi);
146:   VecDuplicate(user.wv,&user.Piv);
147:   VecDuplicate(user.wv,&user.logcv);
148:   VecDuplicate(user.wv,&user.logci);
149:   VecDuplicate(user.wv,&user.logcvi);
150:   VecDuplicate(user.wv,&user.work1);
151:   VecDuplicate(user.wv,&user.work2);
152:   VecDuplicate(user.wv,&user.Riv);
153:   VecDuplicate(user.wv,&user.phi1);
154:   VecDuplicate(user.wv,&user.phi2);
155:   VecDuplicate(user.wv,&user.Phi2D_V);
156:   VecDuplicate(user.wv,&user.Sv);
157:   VecDuplicate(user.wv,&user.Si);
158:   VecDuplicate(user.wv,&user.work3);
159:   VecDuplicate(user.wv,&user.work4);
160:   /* for twodomain modeling */
161:   VecDuplicate(user.wv,&user.phi1);
162:   VecDuplicate(user.wv,&user.phi2);
163:   VecDuplicate(user.wv,&user.Phi2D_V);
164:   VecDuplicate(user.wv,&user.Sv);
165:   VecDuplicate(user.wv,&user.Si);
166:   VecDuplicate(user.wv,&user.work3);
167:   VecDuplicate(user.wv,&user.work4);


170:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
171:   DMCreateMatrix(user.da1,MATAIJ,&user.M);
172:   /* Get the (usual) mass matrix structure from da2 */
173:   DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
174:   SetInitialGuess(x,&user);
175:   /* twodomain modeling */
176:   if (user.domain) {
177:     switch (user.domain) {
178:     case 1:
179:       Phi(&user);
180:       break;
181:     case 2:
182:       Phi_read(&user);
183:       break;
184:     }
185:   }

187:   /* Form the jacobian matrix and M_0 */
188:   SetUpMatrices(&user);
189:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);

191:   SNESCreate(PETSC_COMM_WORLD,&snes);
192:   SNESSetDM(snes,user.da1);

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


198:   SetVariableBounds(user.da1,xl,xu);
199:   SNESVISetVariableBounds(snes,xl,xu);
200:   SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100,PETSC_DEFAULT);
201:   SNESSetFromOptions(snes);

203:   /*  PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_rand",FILE_MODE_WRITE,&view_rand);
204:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat2",FILE_MODE_WRITE,&view_mat);
205:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
206:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
207:    PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);*/

209:   /*PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);

211:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_cv",FILE_MODE_WRITE,&view_cv);
212:    PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_eta",FILE_MODE_WRITE,&view_eta);*/

214:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
215:   if (user.graphics) {
216:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
217:   }
218:   /*
219:   if (user.graphicsfile) {
220:     DMView(user.da1,user.graphicsfile);
221:     VecView(x,user.graphicsfile);
222:   }
223:    */
224:   while (t<user.T) {
225:     SNESSetFunction(snes,r,FormFunction,(void*)&user);
226:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

228:     SetRandomVectors(&user,t);
229:     /*    VecView(user.Pv,view_rand);
230:     VecView(user.Pi,view_rand);
231:      VecView(user.Piv,view_rand);*/

233:     DPsi(&user);
234:     /*    VecView(user.DPsiv,view_psi);
235:     VecView(user.DPsii,view_psi);
236:      VecView(user.DPsieta,view_psi);*/

238:     Update_q(&user);

240:     /*    VecView(user.q,view_q);*/
241:     /*  MatView(user.M,view_mat);*/



245:     sprintf(cv_filename,"file_cv_%f.vtk",t);
246:     sprintf(eta_filename,"file_eta_%f.vtk",t);
247:     /*    PetscViewerASCIIOpen(PETSC_COMM_WORLD,cv_filename,&view_vtk_cv);
248:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,eta_filename,&view_vtk_eta);
249:     PetscViewerSetFormat(view_vtk_cv, PETSC_VIEWER_ASCII_VTK);
250:     PetscViewerSetFormat(view_vtk_eta, PETSC_VIEWER_ASCII_VTK);
251:     DMView(user.da2,view_vtk_cv);
252:     DMView(user.da2,view_vtk_eta);
253:     VecView(user.cv,view_cv);
254:     VecView(user.eta,view_eta);
255:     VecView(user.cv,view_vtk_cv);
256:     VecView(user.eta,view_vtk_eta);
257:     PetscViewerDestroy(&view_vtk_cv);
258:      PetscViewerDestroy(&view_vtk_eta);*/


261:     VecNorm(user.q,NORM_2,&normq);
262:     printf("2-norm of q = %14.12f\n",normq);
263:     SNESSolve(snes,NULL,x);
264:     SNESGetConvergedReason(snes,&reason);
265:     if (reason < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Nonlinear solver failed");
266:     SNESVIGetInactiveSet(snes,&inactiveconstraints);
267:     ISGetSize(inactiveconstraints,&ninactiveconstraints);
268:     /* if (ninactiveconstraints < .90*N) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,"To many active constraints, model has become non-physical"); */

270:     /*    VecView(x,view_out);*/
271:     if (user.graphics) {
272:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
273:     }
274:     /*    VecView(x,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));*/
275:     PetscInt its;
276:     SNESGetIterationNumber(snes,&its);
277:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %g in %d iterations\n",t,its);

279:     Update_u(x,&user);
280:     UpdateMatrices(&user);
281:     t    = t + user.dt;
282:     /*
283:     if (user.graphicsfile) {
284:       VecView(x,user.graphicsfile);
285:     }
286:      */
287:   }

289:   /*  PetscViewerDestroy(&view_rand);
290:   PetscViewerDestroy(&view_mat);
291:   PetscViewerDestroy(&view_q);
292:   PetscViewerDestroy(&view_out);
293:    PetscViewerDestroy(&view_psi);
294:   PetscViewerDestroy(&view_out);

296:   PetscViewerDestroy(&view_cv);
297:    PetscViewerDestroy(&view_eta);*/

299:   if (user.graphicsfile) {
300:     PetscViewerDestroy(&user.graphicsfile);
301:   }

303:   VecDestroy(&x);
304:   VecDestroy(&r);
305:   VecDestroy(&xl);
306:   VecDestroy(&xu);
307:   VecDestroy(&user.q);
308:   VecDestroy(&user.wv);
309:   VecDestroy(&user.cv);
310:   VecDestroy(&user.wi);
311:   VecDestroy(&user.ci);
312:   VecDestroy(&user.eta);
313:   VecDestroy(&user.cvi);
314:   VecDestroy(&user.DPsiv);
315:   VecDestroy(&user.DPsii);
316:   VecDestroy(&user.DPsieta);
317:   VecDestroy(&user.Pv);
318:   VecDestroy(&user.Pi);
319:   VecDestroy(&user.Piv);
320:   VecDestroy(&user.logcv);
321:   VecDestroy(&user.logci);
322:   VecDestroy(&user.logcvi);
323:   VecDestroy(&user.work1);
324:   VecDestroy(&user.work2);
325:   VecDestroy(&user.Riv);
326:   MatDestroy(&user.M);
327:   MatDestroy(&user.M_0);
328:   DMDestroy(&user.da1);
329:   DMDestroy(&user.da2);


332:   PetscFinalize();
333:   return 0;
334: }

338: PetscErrorCode Update_u(Vec X,AppCtx *user)
339: {
341:   PetscInt       i,n;
342:   PetscScalar    *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;

345:   VecGetLocalSize(user->wv,&n);
346:   VecGetArray(X,&xx);
347:   VecGetArray(user->wv,&wv_p);
348:   VecGetArray(user->cv,&cv_p);
349:   VecGetArray(user->wi,&wi_p);
350:   VecGetArray(user->ci,&ci_p);
351:   VecGetArray(user->eta,&eta_p);


354:   for (i=0; i<n; i++) {
355:     wv_p[i]  = xx[5*i];
356:     cv_p[i]  = xx[5*i+1];
357:     wi_p[i]  = xx[5*i+2];
358:     ci_p[i]  = xx[5*i+3];
359:     eta_p[i] = xx[5*i+4];
360:   }
361:   VecRestoreArray(X,&xx);
362:   VecRestoreArray(user->wv,&wv_p);
363:   VecRestoreArray(user->cv,&cv_p);
364:   VecRestoreArray(user->wi,&wi_p);
365:   VecRestoreArray(user->ci,&ci_p);
366:   VecRestoreArray(user->eta,&eta_p);
367:   return(0);
368: }

372: PetscErrorCode Update_q(AppCtx *user)
373: {
375:   PetscScalar    *q_p,*w1,*w2,max1;
376:   PetscInt       i,n;

379:   VecPointwiseMult(user->Riv,user->eta,user->eta);
380:   VecScale(user->Riv,user->Rsurf);
381:   VecShift(user->Riv,user->Rbulk);
382:   VecPointwiseMult(user->Riv,user->ci,user->Riv);
383:   VecPointwiseMult(user->Riv,user->cv,user->Riv);

385:   VecCopy(user->Riv,user->work1);
386:   if (user->radiation) {
387:     VecAXPY(user->work1,-1.0,user->Pv);
388:   }
389:   if (user->domain) {
390:     VecCopy(user->cv,user->work3);
391:     VecShift(user->work3,-1.0*user->cv_eq);
392:     VecCopy(user->Phi2D_V,user->work4);
393:     VecScale(user->work4,-1.0);
394:     VecShift(user->work4,1.0);
395:     VecPointwiseMult(user->work4,user->work4,user->work3);
396:     VecScale(user->work4,user->Svr);
397:     VecAXPY(user->work1,300.0,user->work4);
398:   }
399:   VecScale(user->work1,user->dt);
400:   VecAXPY(user->work1,-1.0,user->cv);
401:   MatMult(user->M_0,user->work1,user->work2);

403:   VecGetArray(user->q,&q_p);
404:   VecGetArray(user->work1,&w1);
405:   VecGetArray(user->work2,&w2);
406:   VecGetLocalSize(user->wv,&n);
407:   for (i=0; i<n; i++) q_p[5*i]=w2[i];

409:   MatMult(user->M_0,user->DPsiv,user->work1);
410:   for (i=0; i<n; i++) q_p[5*i+1]=w1[i];

412:   VecCopy(user->Riv,user->work1);
413:   if (user->radiation) {
414:     VecAXPY(user->work1,-1.0,user->Pi);
415:   }
416:   if (user->domain) {
417:     VecCopy(user->ci,user->work3);
418:     VecShift(user->work3,-1.0*user->ci_eq);
419:     VecCopy(user->Phi2D_V,user->work4);
420:     VecScale(user->work4,-1.0);
421:     VecShift(user->work4,1.0);
422:     VecPointwiseMult(user->work4,user->work4,user->work3);
423:     VecScale(user->work4,user->Sir);
424:     VecAXPY(user->work1,300.0,user->work4);
425:   }
426:   VecScale(user->work1,user->dt);
427:   VecAXPY(user->work1,-1.0,user->ci);
428:   MatMult(user->M_0,user->work1,user->work2);
429:   for (i=0; i<n; i++) q_p[5*i+2]=w2[i];

431:   MatMult(user->M_0,user->DPsii,user->work1);
432:   for (i=0; i<n; i++) q_p[5*i+3]=w1[i];

434:   VecCopy(user->DPsieta,user->work1);
435:   VecScale(user->work1,user->L);
436:   VecAXPY(user->work1,-1.0/user->dt,user->eta);
437:   if (user->radiation) {
438:     VecAXPY(user->work1,-1.0,user->Piv);
439:   }
440:   MatMult(user->M_0,user->work1,user->work2);
441:   for (i=0; i<n; i++) q_p[5*i+4]=w2[i];

443:   VecRestoreArray(user->q,&q_p);
444:   VecRestoreArray(user->work1,&w1);
445:   VecRestoreArray(user->work2,&w2);
446:   return(0);
447: }

451: PetscErrorCode DPsi(AppCtx *user)
452: {
454:   PetscScalar    Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
455:   PetscScalar    *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
456:   PetscInt       n,i;

459:   VecGetLocalSize(user->cv,&n);
460:   VecGetArray(user->cv,&cv_p);
461:   VecGetArray(user->ci,&ci_p);
462:   VecGetArray(user->eta,&eta_p);
463:   VecGetArray(user->logcv,&logcv_p);
464:   VecGetArray(user->logci,&logci_p);
465:   VecGetArray(user->logcvi,&logcvi_p);
466:   VecGetArray(user->DPsiv,&DPsiv_p);
467:   VecGetArray(user->DPsii,&DPsii_p);
468:   VecGetArray(user->DPsieta,&DPsieta_p);

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

472:   Llog(user->ci,user->logci);


475:   VecCopy(user->cv,user->cvi);
476:   VecAXPY(user->cvi,1.0,user->ci);
477:   VecScale(user->cvi,-1.0);
478:   VecShift(user->cvi,1.0);
479:   Llog(user->cvi,user->logcvi);

481:   for (i=0; i<n; i++) {
482:     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);
483:     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];

485:     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]);
486:   }

488:   VecRestoreArray(user->cv,&cv_p);
489:   VecRestoreArray(user->ci,&ci_p);
490:   VecRestoreArray(user->eta,&eta_p);
491:   VecRestoreArray(user->logcv,&logcv_p);
492:   VecRestoreArray(user->logci,&logci_p);
493:   VecRestoreArray(user->logcvi,&logcvi_p);
494:   VecRestoreArray(user->DPsiv,&DPsiv_p);
495:   VecRestoreArray(user->DPsii,&DPsii_p);
496:   VecRestoreArray(user->DPsieta,&DPsieta_p);
497:   return(0);
498: }


503: PetscErrorCode Llog(Vec X, Vec Y)
504: {
506:   PetscScalar    *x,*y;
507:   PetscInt       n,i;

510:   VecGetArray(X,&x);
511:   VecGetArray(Y,&y);
512:   VecGetLocalSize(X,&n);
513:   for (i=0; i<n; i++) {
514:     if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
515:     else y[i] = log(x[i]);
516:   }
517:   return(0);
518: }


523: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
524: {
526:   PetscInt       n,i,Mda,Nda;
527:   PetscScalar    *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta;

529:   /* needed for the void growth case */
530:   PetscScalar       xmid,ymid,cv_v=1.0,cv_m=0.122,ci_v=0.0,ci_m=.00069,eta_v=1.0,eta_m=0.0,h,lambda;
531:   PetscInt          nele,nen,idx[3];
532:   const PetscInt    *ele;
533:   PetscScalar       x[3],y[3];
534:   Vec               coords;
535:   const PetscScalar *_coords;
536:   PetscViewer       view;
537:   PetscScalar       xwidth = user->xmax - user->xmin;

540:   VecGetLocalSize(X,&n);

542:   if (user->voidgrowth) {
543:     DMDAGetInfo(user->da2,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
544:     DMGetCoordinatesLocal(user->da2,&coords);
545:     VecGetArrayRead(coords,&_coords);

547:     h      = (user->xmax-user->xmin)/Mda;
548:     xmid   = (user->xmax + user->xmin)/2.0;
549:     ymid   = (user->ymax + user->ymin)/2.0;
550:     lambda = 4.0*h;

552:     DMDAGetElements(user->da2,&nele,&nen,&ele);
553:     for (i=0; i < nele; i++) {
554:       idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];

556:       x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
557:       x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
558:       x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

560:       PetscInt    k;
561:       PetscScalar vals_cv[3],vals_ci[3],vals_eta[3],s,hhr,r;
562:       for (k=0; k < 3; k++) {
563:         s = PetscSqrtScalar((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
564:         if (s < xwidth*(5.0/64.0)) {
565:           vals_cv[k]  = cv_v;
566:           vals_ci[k]  = ci_v;
567:           vals_eta[k] = eta_v;
568:         } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
569:           /* r = (s - xwidth*(6.0/64.0))/(0.5*lambda); */
570:           r           = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
571:           hhr         = 0.25*(-r*r*r + 3*r + 2);
572:           vals_cv[k]  = cv_m + (1.0 - hhr)*(cv_v - cv_m);
573:           vals_ci[k]  = ci_m + (1.0 - hhr)*(ci_v - ci_m);
574:           vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
575:         } else {
576:           vals_cv[k]  = cv_m;
577:           vals_ci[k]  = ci_m;
578:           vals_eta[k] = eta_m;
579:         }
580:       }
581:       VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
582:       VecSetValuesLocal(user->ci,3,idx,vals_ci,INSERT_VALUES);
583:       VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
584:     }
585:     VecAssemblyBegin(user->cv);
586:     VecAssemblyEnd(user->cv);
587:     VecAssemblyBegin(user->ci);
588:     VecAssemblyEnd(user->ci);
589:     VecAssemblyBegin(user->eta);
590:     VecAssemblyEnd(user->eta);

592:     DMDARestoreElements(user->da2,&nele,&nen,&ele);
593:     VecRestoreArrayRead(coords,&_coords);

595:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view);
596:     VecView(user->cv,view);
597:     VecView(user->ci,view);
598:     VecView(user->eta,view);
599:     PetscViewerDestroy(&view);
600:   } else {
601:     /* VecSet(user->cv,user->initv); */
602:     /* VecSet(user->ci,user->initv); */
603:     VecSet(user->cv,.05);
604:     VecSet(user->ci,.05);
605:     VecSet(user->eta,user->initeta);
606:   }

608:   DPsi(user);
609:   VecCopy(user->DPsiv,user->wv);
610:   VecCopy(user->DPsii,user->wi);

612:   VecGetArray(X,&xx);
613:   VecGetArray(user->cv,&cv_p);
614:   VecGetArray(user->ci,&ci_p);
615:   VecGetArray(user->wv,&wv_p);
616:   VecGetArray(user->wi,&wi_p);
617:   VecGetArray(user->eta,&eta);
618:   for (i=0; i<n/5; i++) {
619:     xx[5*i]  =wv_p[i];
620:     xx[5*i+1]=cv_p[i];
621:     xx[5*i+2]=wi_p[i];
622:     xx[5*i+3]=ci_p[i];
623:     xx[5*i+4]=eta[i];
624:   }

626:   /* VecView(user->wv,view);
627:   VecView(user->cv,view);
628:   VecView(user->wi,view);
629:   VecView(user->ci,view);
630:    PetscViewerDestroy(&view);*/

632:   VecRestoreArray(X,&xx);
633:   VecRestoreArray(user->cv,&cv_p);
634:   VecRestoreArray(user->ci,&ci_p);
635:   VecRestoreArray(user->wv,&wv_p);
636:   VecRestoreArray(user->wi,&wi_p);
637:   VecRestoreArray(user->eta,&eta);
638:   return(0);
639: }

641: typedef struct {
642:   PetscReal dt,x,y,strength;
643: } RandomValues;


648: PetscErrorCode SetRandomVectors(AppCtx *user,PetscReal t)
649: {
650:   PetscErrorCode      ierr;
651:   static RandomValues *randomvalues = 0;
652:   static PetscInt     randindex     = 0,n; /* indicates how far into the randomvalues we have currently used */
653:   static PetscReal    randtime      = 0; /* indicates time of last radiation event */
654:   PetscInt            i,j,M,N,cnt = 0;
655:   PetscInt            xs,ys,xm,ym;

658:   if (!randomvalues) {
659:     PetscViewer viewer;
660:     char        filename[PETSC_MAX_PATH_LEN];
661:     PetscBool   flg;
662:     PetscInt    seed;

664:     PetscOptionsGetInt(NULL,"-random_seed",&seed,&flg);
665:     if (flg) {
666:       sprintf(filename,"ex61.random.%d",(int)seed);
667:     } else {
668:       PetscStrcpy(filename,"ex61.random");
669:     }
670:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
671:     PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
672:     PetscMalloc(n*sizeof(RandomValues),&randomvalues);
673:     PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
674:     for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
675:     PetscViewerDestroy(&viewer);
676:   }
677:   user->maxevents = PetscMin(user->maxevents,n);

679:   VecSet(user->Pv,0.0);
680:   DMDAGetInfo(user->da1,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
681:   DMDAGetGhostCorners(user->da1,&xs,&ys,0,&xm,&ym,0);
682:   while (user->maxevents > randindex && randtime + randomvalues[randindex].dt < t + user->dt) {  /* radiation event has occured since last time step */
683:     i = ((PetscInt) (randomvalues[randindex].x*M)) - xs;
684:     j = ((PetscInt) (randomvalues[randindex].y*N)) - ys;
685:     if (i >= 0 && i < xm && j >= 0 && j < ym) { /* point is on this process */
686:       /* need to make sure eta at the given point is not great than .8 */
687:       VecSetValueLocal(user->Pv,i + 1 + xm*(j + 1), randomvalues[randindex].strength*user->VG,INSERT_VALUES);
688:     }
689:     randtime += randomvalues[randindex++].dt;
690:     cnt++;
691:   }
692:   PetscPrintf(PETSC_COMM_WORLD,"Number of radiation events %d\n",cnt);
693:   VecAssemblyBegin(user->Pv);
694:   VecAssemblyEnd(user->Pv);

696:   VecCopy(user->Pv,user->Pi);
697:   VecScale(user->Pi,0.9);
698:   VecPointwiseMult(user->Piv,user->Pi,user->Pv);
699:   return(0);
700: }

704: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
705: {
707:   AppCtx         *user=(AppCtx*)ctx;

710:   MatMultAdd(user->M,X,user->q,F);
711:   return(0);
712: }

716: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
717: {
719:   AppCtx         *user=(AppCtx*)ctx;

722:   *flg = SAME_NONZERO_PATTERN;
723:   MatCopy(user->M,*J,*flg);
724:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
725:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
726:   return(0);
727: }

731: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
732: {
734:   PetscScalar    ***l,***u;
735:   PetscInt       xs,xm,ys,ym;
736:   PetscInt       j,i;

739:   DMDAVecGetArrayDOF(da,xl,&l);
740:   DMDAVecGetArrayDOF(da,xu,&u);

742:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

744:   for (j=ys; j<ys+ym; j++) {
745:     for (i=xs; i < xs+xm; i++) {
746:       l[j][i][0] = -SNES_VI_INF;
747:       l[j][i][1] = 0.0;
748:       l[j][i][2] = -SNES_VI_INF;
749:       l[j][i][3] = 0.0;
750:       l[j][i][4] = 0.0;
751:       u[j][i][0] = SNES_VI_INF;
752:       u[j][i][1] = 1.0;
753:       u[j][i][2] = SNES_VI_INF;
754:       u[j][i][3] = 1.0;
755:       u[j][i][4] = 1.0;
756:     }
757:   }


760:   DMDAVecRestoreArrayDOF(da,xl,&l);
761:   DMDAVecRestoreArrayDOF(da,xu,&u);
762:   return(0);
763: }

767: PetscErrorCode GetParams(AppCtx *user)
768: {
770:   PetscBool      flg,graphicsfile = PETSC_FALSE;

773:   /* Set default parameters */
774:   user->xmin  = 0.0; user->xmax = 128.0;
775:   user->ymin  = 0.0; user->ymax = 128.0;
776:   user->Dv    = 1.0;
777:   user->Di    = 4.0;
778:   user->Evf   = 0.8;
779:   user->Eif   = 1.2;
780:   user->A     = 1.0;
781:   user->kBT   = 0.11;
782:   user->kav   = 1.0;
783:   user->kai   = 1.0;
784:   user->kaeta = 1.0;
785:   user->Rsurf = 10.0;
786:   user->Rbulk = 1.0;
787:   user->VG    = 100.0;
788:   user->L     = 10.0;

790:   user->T          = 1.0e-2;
791:   user->dt         = 1.0e-4;
792:   user->initv      = .00069;
793:   user->initeta    = 0.0;
794:   user->degenerate = PETSC_FALSE;
795:   user->maxevents  = 10;
796:   user->graphics   = PETSC_TRUE;

798:   /* multidomain modeling */
799:   user->domain =   1;
800:   user->Svr    = 0.5;
801:   user->Sir    = 0.5;
802:   user->cv_eq  = 6.9e-4;
803:   user->ci_eq  = 6.9e-4;
804:   /* void growth */
805:   user->voidgrowth = PETSC_FALSE;

807:   user->radiation = PETSC_FALSE;

809:   /* degenerate mobility */
810:   user->smallnumber = 1.0e-3;

812:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Coupled Cahn-Hillard/Allen-Cahn Equations","Phasefield");
813:   PetscOptionsInt("-domain","Number of domains (0=one domain, 1=two domains, 2=multidomain\n","None",user->domain,&user->domain,&flg);

815:   PetscOptionsReal("-Dv","Vacancy Diffusivity\n","None",user->Dv,&user->Dv,&flg);
816:   PetscOptionsReal("-Di","Interstitial Diffusivity\n","None",user->Di,&user->Di,&flg);
817:   PetscOptionsReal("-Evf","Vacancy Formation Energy\n","None",user->Evf,&user->Evf,&flg);
818:   PetscOptionsReal("-Eif","Interstitial Formation energy\n","None",user->Eif,&user->Eif,&flg);
819:   PetscOptionsReal("-A","???","None",user->A,&user->A,&flg);
820:   PetscOptionsReal("-kBT","Boltzmann's Constant times the Absolute Temperature","None",user->kBT,&user->kBT,&flg);
821:   PetscOptionsReal("-kav","???","None",user->kav,&user->kav,&flg);
822:   PetscOptionsReal("-kai","???","None",user->kai,&user->kai,&flg);
823:   PetscOptionsReal("-kaeta","???","None",user->kaeta,&user->kaeta,&flg);
824:   PetscOptionsReal("-Rsurf","???","None",user->Rsurf,&user->Rsurf,&flg);
825:   PetscOptionsReal("-Rbulk","???","None",user->Rbulk,&user->Rbulk,&flg);
826:   PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
827:   PetscOptionsReal("-L","???","None",user->L,&user->L,&flg);

829:   PetscOptionsReal("-initv","Initial solution of Cv and Ci","None",user->initv,&user->initv,&flg);
830:   PetscOptionsReal("-initeta","Initial solution of Eta","None",user->initeta,&user->initeta,&flg);
831:   PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
832:   PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);

834:   PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
835:   PetscOptionsBool("-radiation","Use initial conditions for void growth\n","None",user->radiation,&user->radiation,&flg);
836:   PetscOptionsReal("-xmin","Lower X coordinate of domain\n","None",user->xmin,&user->xmin,&flg);
837:   PetscOptionsReal("-xmax","Upper X coordinate of domain\n","None",user->xmax,&user->xmax,&flg);
838:   PetscOptionsReal("-T","Total runtime\n","None",user->T,&user->T,&flg);
839:   PetscOptionsReal("-dt","Time step\n","None",user->dt,&user->dt,&flg);

841:   user->dtevent = user->dt;

843:   PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);
844:   PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);

846:   PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
847:   PetscOptionsBool("-graphicsfile","Save solution at each timestep\n","None",graphicsfile,&graphicsfile,&flg);
848:   if (graphicsfile) {
849:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ex61.data",FILE_MODE_WRITE,&user->graphicsfile);
850:   }
851:   PetscOptionsEnd();
852:   return(0);
853: }


858: PetscErrorCode SetUpMatrices(AppCtx *user)
859: {
861:   PetscInt       nele,nen,i,n;
862:   const PetscInt *ele;
863:   PetscScalar    dt=user->dt,hx,hy;

865:   PetscInt    idx[3];
866:   PetscScalar eM_0[3][3],eM_2_even[3][3],eM_2_odd[3][3];
867:   PetscScalar cv_sum, ci_sum;
868:   Mat         M  =user->M;
869:   Mat         M_0=user->M_0;
870:   PetscInt    Mda=user->Mda, Nda=user->Nda;
871:   PetscScalar *cv_p,*ci_p;
872:   /* newly added */
873:   Vec cvlocal,cilocal;

876:   /*  MatSetOption(M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
877:    MatSetOption(M_0,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);*/

879:   /* new stuff */
880:   DMGetLocalVector(user->da2,&cvlocal);
881:   DMGetLocalVector(user->da2,&cilocal);
882:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
883:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
884:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
885:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
886:   /* old stuff */
887:   /*
888:   VecGetArray(user->cv,&cv_p);
889:   VecGetArray(user->ci,&ci_p);
890:    */
891:   /* new stuff */
892:   VecGetArray(cvlocal,&cv_p);
893:   VecGetArray(cilocal,&ci_p);

895:   MatGetLocalSize(M,&n,NULL);
896:   DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
897:   hx   = (user->xmax-user->xmin)/Mda;
898:   hy   = (user->ymax-user->ymin)/Nda;

900:   eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hy/12.0;
901:   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;

903:   eM_2_odd[0][0] = 1.0;
904:   eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
905:   eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
906:   eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;

908:   eM_2_even[0][0] = 1.0;
909:   eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
910:   eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
911:   eM_2_even[1][2] = eM_2_even[2][1] = 0.0;

913:   /*  eM_2_even[1][1] = 1.0;
914:   eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
915:   eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
916:   eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
917:    */

919:   /*  for (k=0;k < Mda*Nda*2;k++) { */
920:   DMDAGetElements(user->da1,&nele,&nen,&ele);
921:   for (i=0; i < nele; i++) {
922:     /*
923:     idx[0] = connect[k*3];
924:     idx[1] = connect[k*3+1];
925:     idx[2] = connect[k*3+2];
926:      */
927:     idx[0] = ele[3*i];
928:     idx[1] = ele[3*i+1];
929:     idx[2] = ele[3*i+2];

931:     PetscInt    row,cols[6],r,row_M_0,cols3[3];
932:     PetscScalar vals[6],vals_M_0[3],vals3[3];

934:     for (r=0; r<3; r++) {
935:       /* row_M_0 = connect[k*3+r]; */
936:       row_M_0 = idx[r];

938:       vals_M_0[0]=eM_0[r][0];
939:       vals_M_0[1]=eM_0[r][1];
940:       vals_M_0[2]=eM_0[r][2];


943:       MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);

945:       if (user->degenerate) {
946:         cv_sum = (cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
947:         ci_sum = (ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
948:       } else {
949:         cv_sum = user->initv*user->Dv/user->kBT;
950:         ci_sum = user->initv*user->Di/user->kBT;
951:       }



955:       row     = 5*idx[r];
956:       cols[0] = 5*idx[0];     vals[0] = dt*eM_2_odd[r][0]*cv_sum;
957:       cols[1] = 5*idx[1];     vals[1] = dt*eM_2_odd[r][1]*cv_sum;
958:       cols[2] = 5*idx[2];     vals[2] = dt*eM_2_odd[r][2]*cv_sum;
959:       cols[3] = 5*idx[0]+1;   vals[3] = eM_0[r][0];
960:       cols[4] = 5*idx[1]+1;   vals[4] = eM_0[r][1];
961:       cols[5] = 5*idx[2]+1;   vals[5] = eM_0[r][2];


964:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);


967:       row     = 5*idx[r]+1;
968:       cols[0] = 5*idx[0];     vals[0] = -1.0*eM_0[r][0];
969:       cols[1] = 5*idx[1];     vals[1] = -1.0*eM_0[r][1];
970:       cols[2] = 5*idx[2];     vals[2] = -1.0*eM_0[r][2];
971:       cols[3] = 5*idx[0]+1;   vals[3] =  user->kav*eM_2_odd[r][0];
972:       cols[4] = 5*idx[1]+1;   vals[4] =  user->kav*eM_2_odd[r][1];
973:       cols[5] = 5*idx[2]+1;   vals[5] =  user->kav*eM_2_odd[r][2];

975:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);

977:       row     = 5*idx[r]+2;
978:       cols[0] = 5*idx[0]+2;   vals[0] =  dt*eM_2_odd[r][0]*ci_sum;
979:       cols[1] = 5*idx[1]+2;   vals[1] =  dt*eM_2_odd[r][1]*ci_sum;
980:       cols[2] = 5*idx[2]+2;   vals[2] =  dt*eM_2_odd[r][2]*ci_sum;
981:       cols[3] = 5*idx[0]+3;   vals[3] =  eM_0[r][0];
982:       cols[4] = 5*idx[1]+3;   vals[4] =  eM_0[r][1];
983:       cols[5] = 5*idx[2]+3;   vals[5] =  eM_0[r][2];

985:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);


988:       row     = 5*idx[r]+3;
989:       cols[0] = 5*idx[0]+2;   vals[0] = -1.0*eM_0[r][0];
990:       cols[1] = 5*idx[1]+2;   vals[1] = -1.0*eM_0[r][1];
991:       cols[2] = 5*idx[2]+2;   vals[2] = -1.0*eM_0[r][2];
992:       cols[3] = 5*idx[0]+3;   vals[3] =  user->kai*eM_2_odd[r][0];
993:       cols[4] = 5*idx[1]+3;   vals[4] =  user->kai*eM_2_odd[r][1];
994:       cols[5] = 5*idx[2]+3;   vals[5] =  user->kai*eM_2_odd[r][2];

996:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);


999:       row = 5*idx[r]+4;
1000:       /*
1001:       cols3[0] = 5*idx[0]+4;   vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_odd[r][0];
1002:       cols3[1] = 5*idx[1]+4;   vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_odd[r][1];
1003:       cols3[2] = 5*idx[2]+4;   vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_odd[r][2];
1004:        */
1005:       cols3[0] = 5*idx[0]+4;   vals3[0] = (eM_0[r][0]/dt + user->L*user->kaeta*eM_2_odd[r][0]);
1006:       cols3[1] = 5*idx[1]+4;   vals3[1] = (eM_0[r][1]/dt + user->L*user->kaeta*eM_2_odd[r][1]);
1007:       cols3[2] = 5*idx[2]+4;   vals3[2] = (eM_0[r][2]/dt + user->L*user->kaeta*eM_2_odd[r][2]);

1009:       MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);

1011:     }
1012:   }

1014:   /* new */
1015:   VecRestoreArray(cvlocal,&cv_p);
1016:   VecRestoreArray(cilocal,&ci_p);
1017:   DMRestoreLocalVector(user->da2,&cvlocal);
1018:   DMRestoreLocalVector(user->da2,&cilocal);
1019:   /* old */
1020:   /*
1021:   VecRestoreArray(user->cv,&cv_p);
1022:   VecRestoreArray(user->ci,&ci_p);
1023:    */
1024:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
1025:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);

1027:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1028:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

1030:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
1031:   return(0);
1032: }


1037: PetscErrorCode UpdateMatrices(AppCtx *user)
1038: {
1040:   PetscInt       i,n,Mda,Nda,nele,nen;
1041:   const PetscInt *ele;

1043:   PetscInt    idx[3];
1044:   PetscScalar eM_2_odd[3][3],eM_2_even[3][3],h,dt=user->dt;
1045:   Mat         M=user->M;
1046:   PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
1047:   /* newly added */
1048:   Vec cvlocal,cilocal;

1051:   MatGetLocalSize(M,&n,NULL);

1053:   /* new stuff */
1054:   DMGetLocalVector(user->da2,&cvlocal);
1055:   DMGetLocalVector(user->da2,&cilocal);
1056:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
1057:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
1058:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
1059:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
1060:   /* new stuff */
1061:   VecGetArray(cvlocal,&cv_p);
1062:   VecGetArray(cilocal,&ci_p);

1064:   /* old stuff */
1065:   /*
1066:   VecGetArray(user->cv,&cv_p);
1067:   VecGetArray(user->ci,&ci_p);
1068:    */
1069:   DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);



1073:   h = (user->xmax-user->xmin)/Mda;

1075:   DMDAGetElements(user->da1,&nele,&nen,&ele);

1077:   for (i=0; i < nele; i++) {
1078:     /*
1079:     idx[0] = connect[k*3];
1080:     idx[1] = connect[k*3+1];
1081:     idx[2] = connect[k*3+2];
1082:      */
1083:     idx[0] = ele[3*i];
1084:     idx[1] = ele[3*i+1];
1085:     idx[2] = ele[3*i+2];

1087:     PetscInt    r,row,cols[3];
1088:     PetscScalar vals[3];
1089:     for (r=0; r<3; r++) {
1090:       row     = 5*idx[r];
1091:       cols[0] = 5*idx[0];     vals[0] = 0.0;
1092:       cols[1] = 5*idx[1];     vals[1] = 0.0;
1093:       cols[2] = 5*idx[2];     vals[2] = 0.0;

1095:       /* Insert values in matrix M for 1st dof */
1096:       MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);

1098:       row     = 5*idx[r]+2;
1099:       cols[0] = 5*idx[0]+2;   vals[0] = 0.0;
1100:       cols[1] = 5*idx[1]+2;   vals[1] = 0.0;
1101:       cols[2] = 5*idx[2]+2;   vals[2] = 0.0;

1103:       /* Insert values in matrix M for 3nd dof */
1104:       MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1105:     }
1106:   }

1108:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1109:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

1111:   eM_2_odd[0][0] = 1.0;
1112:   eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
1113:   eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
1114:   eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;

1116:   eM_2_even[0][0] = 1.0;
1117:   eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
1118:   eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
1119:   eM_2_even[1][2] = eM_2_even[2][1] = 0.0;

1121:   /*
1122:   eM_2_even[1][1] = 1.0;
1123:   eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
1124:   eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
1125:   eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
1126:    */

1128:   /* Get local element info */
1129:   /* for (k=0;k < Mda*Nda*2;k++) { */
1130:   for (i=0; i < nele; i++) {
1131:     /*
1132:       idx[0] = connect[k*3];
1133:       idx[1] = connect[k*3+1];
1134:       idx[2] = connect[k*3+2];
1135:      */
1136:     idx[0] = ele[3*i];
1137:     idx[1] = ele[3*i+1];
1138:     idx[2] = ele[3*i+2];

1140:     PetscInt    row,cols[3],r;
1141:     PetscScalar vals[3];

1143:     for (r=0; r<3; r++) {

1145:       if (user->degenerate) {
1146:         printf("smallnumber = %14.12f\n",user->smallnumber);
1147:         cv_sum = (user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
1148:         ci_sum = (user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
1149:       } else {
1150:         cv_sum = user->initv*user->Dv/(user->kBT);
1151:         ci_sum = user->initv*user->Di/user->kBT;
1152:       }

1154:       row     = 5*idx[r];
1155:       cols[0] = 5*idx[0];     vals[0] = dt*eM_2_odd[r][0]*cv_sum;
1156:       cols[1] = 5*idx[1];     vals[1] = dt*eM_2_odd[r][1]*cv_sum;
1157:       cols[2] = 5*idx[2];     vals[2] = dt*eM_2_odd[r][2]*cv_sum;

1159:       /* Insert values in matrix M for 1st dof */
1160:       MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);

1162:       row     = 5*idx[r]+2;
1163:       cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2_odd[r][0]*ci_sum;
1164:       cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2_odd[r][1]*ci_sum;
1165:       cols[2] = 5*idx[2]+2;   vals[2] = dt*eM_2_odd[r][2]*ci_sum;

1167:       MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
1168:     }
1169:   }

1171:   /* new stuff */
1172:   VecRestoreArray(cvlocal,&cv_p);
1173:   VecRestoreArray(cilocal,&ci_p);
1174:   DMRestoreLocalVector(user->da2,&cvlocal);
1175:   DMRestoreLocalVector(user->da2,&cilocal);
1176:   /* old stuff */
1177:   /*
1178:   VecRestoreArray(user->cv,&cv_p);
1179:   VecRestoreArray(user->ci,&ci_p);
1180:    */

1182:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1183:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1184:   return(0);
1185: }



1191: PetscErrorCode Phi(AppCtx *user)
1192: {
1193:   PetscErrorCode    ierr;
1194:   PetscScalar       xmid, xqu, lambda, h,x[3],y[3];
1195:   Vec               coords;
1196:   const PetscScalar *_coords;
1197:   PetscInt          nele,nen,i,idx[3],Mda,Nda;
1198:   const PetscInt    *ele;
1199:   PetscViewer       view;

1202:   DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1203:   DMGetCoordinatesLocal(user->da2,&coords);
1204:   VecGetArrayRead(coords,&_coords);

1206:   h      = (user->xmax - user->xmin)/Mda;
1207:   xmid   = (user->xmin + user->xmax)/2.0;
1208:   xqu    = (user->xmin + xmid)/2.0;
1209:   lambda = 4.0*h;


1212:   DMDAGetElements(user->da2,&nele,&nen,&ele);
1213:   for (i=0; i < nele; i++) {
1214:     idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
1215:     /* printf("idx[0]=%d,idx[1]=%d,idx[2]=%d\n",idx[0],idx[1],idx[2]); */

1217:     x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
1218:     x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
1219:     x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

1221:     /*printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);*/
1222:     /*printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);*/

1224:     PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
1225:     PetscInt    k;

1227:     xc1 = user->xmin;
1228:     xc2 = xmid;

1230:     /*VecSet(user->phi1,0.0);*/
1231:     for (k=0; k < 3; k++) {
1232:       if (x[k]-xqu > 0) s1 = (x[k] - xqu);
1233:       else s1 = -(x[k] - xqu);

1235:       if (x[k] - xc1 > 0) dist1 = (x[k] - xc1);
1236:       else dist1 = -(x[k] - xc1);

1238:       if (x[k] - xc2 > 0) dist2 = (x[k] - xc2);
1239:       else dist2 = -(x[k] - xc2);

1241:       if (dist1 <= 0.5*lambda) {
1242:         r        = (x[k]-xc1)/(0.5*lambda);
1243:         hhr      = 0.25*(-r*r*r + 3.0*r + 2.0);
1244:         vals1[k] = hhr;
1245:       } else if (dist2 <= 0.5*lambda) {
1246:         r        = (x[k]-xc2)/(0.5*lambda);
1247:         hhr      = 0.25*(-r*r*r + 3.0*r + 2.0);
1248:         vals1[k] = 1.0 - hhr;
1249:       } else if (s1 <= xqu - 2.0*h) vals1[k] = 1.0;
1250:       /*else if (abs(x[k]-(user->xmax-h)) < 0.1*h) {*/
1251:       else if ((user->xmax-h)-x[k] < 0.1*h) vals1[k] = .15625;
1252:       else vals1[k] = 0.0;
1253:     }

1255:     VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);

1257:     xc1 = xmid;
1258:     xc2 = user->xmax;

1260:     /*VecSet(user->phi2,0.0);*/
1261:     for (k=0;k < 3; k++) {
1262:       /*
1263:       s1 = abs(x[k] - (xqu+xmid));
1264:       dist1 = abs(x[k] - xc1);
1265:       dist2 = abs(x[k] - xc2);
1266:        */
1267:       if (x[k]-(xqu + xmid) > 0) s1 = (x[k] - (xqu + xmid));
1268:       else s1 = -(x[k] - (xqu + xmid));

1270:       if (x[k] - xc1 > 0) dist1 = (x[k] - xc1);
1271:       else dist1 = -(x[k] - xc1);

1273:       if (x[k] - xc2 > 0) dist2 = (x[k] - xc2);
1274:       else dist2 = -(x[k] - xc2);

1276:       if (dist1 <= 0.5*lambda) {
1277:         r        = (x[k]-xc1)/(0.5*lambda);
1278:         hhr      = 0.25*(-r*r*r + 3.0*r + 2.0);
1279:         vals2[k] = hhr;
1280:       } else if (dist2 <= 0.5*lambda) {
1281:         r        = -(x[k]-xc2)/(0.5*lambda);
1282:         hhr      = 0.25*(-r*r*r + 3.0*r + 2.0);
1283:         vals2[k] = hhr;
1284:       } else if (s1 <= xqu - 2.0*h)           vals2[k] = 1.0;
1285:       else if (x[k]-(user->xmin) < 0.1*h)     vals2[k] = 0.5;
1286:       else if ((x[k]-(user->xmin+h)) < 0.1*h) vals2[k] = .15625;
1287:       else vals2[k] = 0.0;

1289:     }

1291:     VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
1292:     /*
1293:     for (k=0;k < 3; k++) {
1294:       vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
1295:     }
1296:      */
1297:     /*VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);*/

1299:   }

1301:   VecAssemblyBegin(user->phi1);
1302:   VecAssemblyEnd(user->phi1);
1303:   VecAssemblyBegin(user->phi2);
1304:   VecAssemblyEnd(user->phi2);
1305:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);

1307:   VecView(user->phi1,view);
1308:   VecView(user->phi2,view);


1311:   /*VecView(user->phi1,0);*/
1312:   /*VecView(user->phi2,0);*/

1314:   VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1315:   VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1316:   VecView(user->phi1,view);
1317:   VecView(user->phi2,view);

1319:   VecCopy(user->phi1,user->Phi2D_V);
1320:   VecAXPY(user->Phi2D_V,1.0,user->phi2);
1321:   /*VecView(user->Phi2D_V,0);*/

1323:   VecView(user->Phi2D_V,view);
1324:   PetscViewerDestroy(&view);
1325:   /*  VecNorm(user->Phi2D_V,NORM_INFINITY,&max1);*/
1326:   /*VecMin(user->Phi2D_V,&loc1,&min1);*/
1327:   /*printf("norm phi = %f, min phi = %f\n",max1,min1);*/
1328:   return(0);
1329: }

1333: PetscErrorCode Phi_read(AppCtx *user)
1334: {
1336:   PetscReal      *values;
1337:   PetscViewer    viewer;

1340:   VecGetArray(user->Phi2D_V,&values);
1341:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_READ,&viewer);
1342:   PetscViewerBinaryRead(viewer,values,16384,PETSC_DOUBLE);
1343:   PetscViewerDestroy(&viewer);
1344:   return(0);
1345: }