Actual source code: ex61.c

  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_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_ls 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_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_ls 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_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_ls 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_ls_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_ls_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_ls 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

 38:    
 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:   PetscBool   twodomain;
 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*);

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

105:   PetscInitialize(&argc,&argv, (char*)0, help);
106: 
107:   /* Get physics and time parameters */
108:   GetParams(&user);
109:   /* Create a 1D DA with dof = 5; the whole thing */
110:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 5, 1,PETSC_NULL,PETSC_NULL,&user.da1);
111: 
112:   /* Create a 1D DA with dof = 1; for individual componentes */
113:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 1, 1,PETSC_NULL,PETSC_NULL,&user.da2);


116:   /* Set Element type (triangular) */
117:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
118:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
119: 
120:   /* Set x and y coordinates */
121:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
122:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);


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


169:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
170:   DMGetMatrix(user.da1,MATAIJ,&user.M);
171:   /* Get the (usual) mass matrix structure from da2 */
172:   DMGetMatrix(user.da2,MATAIJ,&user.M_0);
173:   SetInitialGuess(x,&user);
174:   /* twodomain modeling */
175:   if (user.twodomain) {
176:     Phi(&user);
177:   }

179:   /* Form the jacobian matrix and M_0 */
180:   SetUpMatrices(&user);
181:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
182: 
183:   SNESCreate(PETSC_COMM_WORLD,&snes);
184:   SNESSetDM(snes,user.da1);

186:   SNESSetFunction(snes,r,FormFunction,(void*)&user);
187:   SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
188: 

190:   SNESSetType(snes,SNESVI);
191:   SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100,PETSC_DEFAULT);

193:   SNESSetFromOptions(snes);
194:   SetVariableBounds(user.da1,xl,xu);
195:   SNESVISetVariableBounds(snes,xl,xu);
196: 
197:   /*  PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_rand",FILE_MODE_WRITE,&view_rand);
198:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat2",FILE_MODE_WRITE,&view_mat);
199:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
200:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
201:    PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);*/
202: 
203:   /*PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
204:   
205:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_cv",FILE_MODE_WRITE,&view_cv);
206:    PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_eta",FILE_MODE_WRITE,&view_eta);*/
207: 
208:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
209:   if (user.graphics) {
210:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
211:   }
212:   if (user.graphicsfile) {
213:     DMView(user.da1,user.graphicsfile);
214:     VecView(x,user.graphicsfile);
215:   }
216:   while (t<user.T) {
217:     SNESSetFunction(snes,r,FormFunction,(void*)&user);
218:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

220:     SetRandomVectors(&user,t);
221:     /*    VecView(user.Pv,view_rand);
222:     VecView(user.Pi,view_rand);
223:      VecView(user.Piv,view_rand);*/

225:     DPsi(&user);
226:     /*    VecView(user.DPsiv,view_psi);
227:     VecView(user.DPsii,view_psi);
228:      VecView(user.DPsieta,view_psi);*/

230:     Update_q(&user);

232:     /*    VecView(user.q,view_q);*/
233:     /*  MatView(user.M,view_mat);*/

235: 
236: 
237:     sprintf(cv_filename,"file_cv_%f.vtk",t);
238:     sprintf(eta_filename,"file_eta_%f.vtk",t);
239:     /*    PetscViewerASCIIOpen(PETSC_COMM_WORLD,cv_filename,&view_vtk_cv);
240:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,eta_filename,&view_vtk_eta);
241:     PetscViewerSetFormat(view_vtk_cv, PETSC_VIEWER_ASCII_VTK);
242:     PetscViewerSetFormat(view_vtk_eta, PETSC_VIEWER_ASCII_VTK);
243:     DMView(user.da2,view_vtk_cv);
244:     DMView(user.da2,view_vtk_eta);
245:     VecView(user.cv,view_cv);
246:     VecView(user.eta,view_eta);
247:     VecView(user.cv,view_vtk_cv);
248:     VecView(user.eta,view_vtk_eta);
249:     PetscViewerDestroy(&view_vtk_cv);
250:      PetscViewerDestroy(&view_vtk_eta);*/

252: 
253:     VecNorm(user.q,NORM_2,&normq);
254:     printf("2-norm of q = %14.12f\n",normq);
255:     SNESSolve(snes,PETSC_NULL,x);
256:     SNESGetConvergedReason(snes,&reason);
257:     if (reason < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Nonlinear solver failed");
258:     SNESVIGetInactiveSet(snes,&inactiveconstraints);
259:     ISGetSize(inactiveconstraints,&ninactiveconstraints);
260:     /* if (ninactiveconstraints < .90*N) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,"To many active constraints, model has become non-physical"); */

262:     /*    VecView(x,view_out);*/
263:     if (user.graphics) {
264:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
265:     }
266:     /*    VecView(x,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));*/
267:     PetscInt its;
268:     SNESGetIterationNumber(snes,&its);
269:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %g in %d iterations\n",t,its);

271:     Update_u(x,&user);
272:     UpdateMatrices(&user);
273:     t = t + user.dt;
274:     if (user.graphicsfile) {
275:       VecView(x,user.graphicsfile);
276:     }
277:   }
278: 
279:   /*  PetscViewerDestroy(&view_rand);
280:   PetscViewerDestroy(&view_mat);
281:   PetscViewerDestroy(&view_q);
282:   PetscViewerDestroy(&view_out);
283:    PetscViewerDestroy(&view_psi);
284:   PetscViewerDestroy(&view_out);
285:   
286:   PetscViewerDestroy(&view_cv);
287:    PetscViewerDestroy(&view_eta);*/
288: 
289:   if (user.graphicsfile) {
290:     PetscViewerDestroy(&user.graphicsfile);
291:   }
292: 
293:   VecDestroy(&x);
294:   VecDestroy(&r);
295:   VecDestroy(&xl);
296:   VecDestroy(&xu);
297:   VecDestroy(&user.q);
298:   VecDestroy(&user.wv);
299:   VecDestroy(&user.cv);
300:   VecDestroy(&user.wi);
301:   VecDestroy(&user.ci);
302:   VecDestroy(&user.eta);
303:   VecDestroy(&user.cvi);
304:   VecDestroy(&user.DPsiv);
305:   VecDestroy(&user.DPsii);
306:   VecDestroy(&user.DPsieta);
307:   VecDestroy(&user.Pv);
308:   VecDestroy(&user.Pi);
309:   VecDestroy(&user.Piv);
310:   VecDestroy(&user.logcv);
311:   VecDestroy(&user.logci);
312:   VecDestroy(&user.logcvi);
313:   VecDestroy(&user.work1);
314:   VecDestroy(&user.work2);
315:   VecDestroy(&user.Riv);
316:   MatDestroy(&user.M);
317:   MatDestroy(&user.M_0);
318:   DMDestroy(&user.da1);
319:   DMDestroy(&user.da2);

321: 
322:   PetscFinalize();
323:   return 0;
324: }

328: PetscErrorCode Update_u(Vec X,AppCtx *user)
329: {
331:   PetscInt       i,n;
332:   PetscScalar    *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
333: 
335:   VecGetLocalSize(user->wv,&n);
336:   VecGetArray(X,&xx);
337:   VecGetArray(user->wv,&wv_p);
338:   VecGetArray(user->cv,&cv_p);
339:   VecGetArray(user->wi,&wi_p);
340:   VecGetArray(user->ci,&ci_p);
341:   VecGetArray(user->eta,&eta_p);


344:   for(i=0;i<n;i++) {
345:     wv_p[i] = xx[5*i];
346:     cv_p[i] = xx[5*i+1];
347:     wi_p[i] = xx[5*i+2];
348:     ci_p[i] = xx[5*i+3];
349:     eta_p[i] = xx[5*i+4];
350:   }
351:   VecRestoreArray(X,&xx);
352:   VecRestoreArray(user->wv,&wv_p);
353:   VecRestoreArray(user->cv,&cv_p);
354:   VecRestoreArray(user->wi,&wi_p);
355:   VecRestoreArray(user->ci,&ci_p);
356:   VecRestoreArray(user->eta,&eta_p);
357: 
358:   return(0);
359: }

363: PetscErrorCode Update_q(AppCtx *user)
364: {
366:   PetscScalar    *q_p,*w1,*w2,max1;
367:   PetscInt       i,n;

369: 
371: 
372:   VecPointwiseMult(user->Riv,user->eta,user->eta);
373:   VecScale(user->Riv,user->Rsurf);
374:   VecShift(user->Riv,user->Rbulk);
375:   VecPointwiseMult(user->Riv,user->ci,user->Riv);
376:   VecPointwiseMult(user->Riv,user->cv,user->Riv);
377: 
378:   VecCopy(user->Riv,user->work1);
379:   VecScale(user->work1,user->dt);
380:   VecAXPY(user->work1,-1.0,user->cv);
381:   MatMult(user->M_0,user->work1,user->work2);
382: 
383:   VecGetArray(user->q,&q_p);
384:   VecGetArray(user->work1,&w1);
385:   VecGetArray(user->work2,&w2);
386:   VecGetLocalSize(user->wv,&n);
387:   for (i=0;i<n;i++) {
388:     q_p[5*i]=w2[i];
389:   }

391:   MatMult(user->M_0,user->DPsiv,user->work1);
392:   for (i=0;i<n;i++) {
393:     q_p[5*i+1]=w1[i];
394:   }

396:   VecCopy(user->Riv,user->work1);
397:   VecScale(user->work1,user->dt);
398:   VecAXPY(user->work1,-1.0,user->ci);
399:   MatMult(user->M_0,user->work1,user->work2);
400:   for (i=0;i<n;i++) {
401:     q_p[5*i+2]=w2[i];
402:   }

404:   MatMult(user->M_0,user->DPsii,user->work1);
405:   for (i=0;i<n;i++) {
406:     q_p[5*i+3]=w1[i];
407:   }

409:   VecCopy(user->DPsieta,user->work1);
410:   VecScale(user->work1,user->L);
411:   VecAXPY(user->work1,-1.0/user->dt,user->eta);
412:   MatMult(user->M_0,user->work1,user->work2);
413:   for (i=0;i<n;i++) {
414:     q_p[5*i+4]=w2[i];
415:   }
416: 
417:   VecRestoreArray(user->q,&q_p);
418:   VecRestoreArray(user->work1,&w1);
419:   VecRestoreArray(user->work2,&w2);
420: 
421:   return(0);
422: }

426: PetscErrorCode DPsi(AppCtx* user)
427: {
428:   PetscErrorCode  ierr;
429:   PetscScalar     Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
430:   PetscScalar     *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
431:   PetscInt        n,i;


435:   VecGetLocalSize(user->cv,&n);
436:   VecGetArray(user->cv,&cv_p);
437:   VecGetArray(user->ci,&ci_p);
438:   VecGetArray(user->eta,&eta_p);
439:   VecGetArray(user->logcv,&logcv_p);
440:   VecGetArray(user->logci,&logci_p);
441:   VecGetArray(user->logcvi,&logcvi_p);
442:   VecGetArray(user->DPsiv,&DPsiv_p);
443:   VecGetArray(user->DPsii,&DPsii_p);
444:   VecGetArray(user->DPsieta,&DPsieta_p);
445: 
446:   Llog(user->cv,user->logcv);
447: 
448:   Llog(user->ci,user->logci);
449: 

451:   VecCopy(user->cv,user->cvi);
452:   VecAXPY(user->cvi,1.0,user->ci);
453:   VecScale(user->cvi,-1.0);
454:   VecShift(user->cvi,1.0);
455:   Llog(user->cvi,user->logcvi);
456: 
457:   for (i=0;i<n;i++)
458:   {
459:     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);

461:     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] ;
462: 
463:     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]);
464: 
465: 
466:   }

468:   VecRestoreArray(user->cv,&cv_p);
469:   VecRestoreArray(user->ci,&ci_p);
470:   VecRestoreArray(user->eta,&eta_p);
471:   VecRestoreArray(user->logcv,&logcv_p);
472:   VecRestoreArray(user->logci,&logci_p);
473:   VecRestoreArray(user->logcvi,&logcvi_p);
474:   VecRestoreArray(user->DPsiv,&DPsiv_p);
475:   VecRestoreArray(user->DPsii,&DPsii_p);
476:   VecRestoreArray(user->DPsieta,&DPsieta_p);


479:   return(0);

481: }


486: PetscErrorCode Llog(Vec X, Vec Y)
487: {
488:   PetscErrorCode    ierr;
489:   PetscScalar       *x,*y;
490:   PetscInt          n,i;

493: 
494:   VecGetArray(X,&x);
495:   VecGetArray(Y,&y);
496:   VecGetLocalSize(X,&n);
497:   for (i=0;i<n;i++) {
498:     if (x[i] < 1.0e-12) {
499:       y[i] = log(1.0e-12);
500:     }
501:     else {
502:       y[i] = log(x[i]);
503:     }
504:   }
505: 
506:   return(0);
507: }


512: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
513: {
514:   PetscErrorCode    ierr;
515:   PetscInt          n,i,Mda,Nda;
516:   PetscScalar           *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta;

518:   /* needed for the void growth case */
519:   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;
520:   PetscInt          nele,nen,idx[3];
521:   const PetscInt    *ele;
522:   PetscScalar       x[3],y[3];
523:   Vec               coords;
524:   const PetscScalar *_coords;
525:   PetscViewer       view;
526:   PetscScalar       xwidth = user->xmax - user->xmin;


530:   VecGetLocalSize(X,&n);

532:   if (user->voidgrowth) {
533:     DMDAGetInfo(user->da2,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
534:     DMDAGetGhostedCoordinates(user->da2,&coords);
535:     VecGetArrayRead(coords,&_coords);

537:     h = (user->xmax-user->xmin)/Mda;
538:     xmid = (user->xmax + user->xmin)/2.0;
539:     ymid = (user->ymax + user->ymin)/2.0;
540:     lambda = 4.0*h;

542:     DMDAGetElements(user->da2,&nele,&nen,&ele);
543:     for (i=0;i < nele; i++) {
544:       idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
545: 
546:       x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
547:       x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
548:       x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
549: 
550:       PetscInt k;
551:       PetscScalar vals_cv[3],vals_ci[3],vals_eta[3],s,hhr,r;
552:       for (k=0; k < 3 ; k++) {
553:         s = sqrt((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
554:         if (s < xwidth*(5.0/64.0)) {
555:           vals_cv[k] = cv_v;
556:           vals_ci[k] = ci_v;
557:           vals_eta[k] = eta_v;
558:         } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
559:           //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
560:           r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
561:           hhr = 0.25*(-r*r*r + 3*r + 2);
562:           vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
563:           vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
564:           vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
565:         } else {
566:           vals_cv[k] = cv_m;
567:           vals_ci[k] = ci_m;
568:           vals_eta[k] = eta_m;
569:         }
570:       }
571:       VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
572:       VecSetValuesLocal(user->ci,3,idx,vals_ci,INSERT_VALUES);
573:       VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
574:     }
575:     VecAssemblyBegin(user->cv);
576:     VecAssemblyEnd(user->cv);
577:     VecAssemblyBegin(user->ci);
578:     VecAssemblyEnd(user->ci);
579:     VecAssemblyBegin(user->eta);
580:     VecAssemblyEnd(user->eta);

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

585:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view);
586:     VecView(user->cv,view);
587:     VecView(user->ci,view);
588:     VecView(user->eta,view);
589:     PetscViewerDestroy(&view);
590:   }
591:   else {
592:     VecSet(user->cv,user->initv);
593:     VecSet(user->ci,user->initv);
594:     VecSet(user->eta,user->initeta);
595:   }

597:   DPsi(user);
598:   VecCopy(user->DPsiv,user->wv);
599:   VecCopy(user->DPsii,user->wi);

601:   VecGetArray(X,&xx);
602:   VecGetArray(user->cv,&cv_p);
603:   VecGetArray(user->ci,&ci_p);
604:   VecGetArray(user->wv,&wv_p);
605:   VecGetArray(user->wi,&wi_p);
606:   VecGetArray(user->eta,&eta);
607:   for (i=0;i<n/5;i++)
608:   {
609:     xx[5*i]=wv_p[i];
610:     xx[5*i+1]=cv_p[i];
611:     xx[5*i+2]=wi_p[i];
612:     xx[5*i+3]=ci_p[i];
613:     xx[5*i+4]=eta[i];
614:   }

616:   /* VecView(user->wv,view);
617:   VecView(user->cv,view);
618:   VecView(user->wi,view);
619:   VecView(user->ci,view);
620:    PetscViewerDestroy(&view);*/

622:   VecRestoreArray(X,&xx);
623:   VecRestoreArray(user->cv,&cv_p);
624:   VecRestoreArray(user->ci,&ci_p);
625:   VecRestoreArray(user->wv,&wv_p);
626:   VecRestoreArray(user->wi,&wi_p);
627:   VecRestoreArray(user->eta,&eta);
628: 
629:   return(0);
630: }

632: typedef struct {
633:   PetscReal dt,x,y,strength;
634: } RandomValues;


639: PetscErrorCode SetRandomVectors(AppCtx* user,PetscReal t)
640: {
641:   PetscErrorCode        ierr;
642:   static RandomValues   *randomvalues = 0;
643:   static PetscInt       randindex = 0,n; /* indicates how far into the randomvalues we have currently used */
644:   static PetscReal      randtime = 0; /* indicates time of last radiation event */
645:   PetscInt              i,j,M,N,cnt = 0;
646:   PetscInt              xs,ys,xm,ym;

649:   if (!randomvalues) {
650:     PetscViewer viewer;
651:     char        filename[PETSC_MAX_PATH_LEN];
652:     PetscBool   flg;
653:     PetscInt    seed;

655:     PetscOptionsGetInt(PETSC_NULL,"-random_seed",&seed,&flg);
656:     if (flg) {
657:       sprintf(filename,"ex61.random.%d",(int)seed);
658:     } else {
659:       PetscStrcpy(filename,"ex61.random");
660:     }
661:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
662:     PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
663:     PetscMalloc(n*sizeof(RandomValues),&randomvalues);
664:     PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
665:     for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
666:     PetscViewerDestroy(&viewer);
667:   }
668:   user->maxevents = PetscMin(user->maxevents,n);

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

687:   VecCopy(user->Pv,user->Pi);
688:   VecScale(user->Pi,0.9);
689:   VecPointwiseMult(user->Piv,user->Pi,user->Pv);
690:   return(0);
691: 
692: }

696: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
697: {
699:   AppCtx         *user=(AppCtx*)ctx;
700: 
702:   MatMultAdd(user->M,X,user->q,F);
703:   return(0);
704: }

708: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
709: {
711:   AppCtx         *user=(AppCtx*)ctx;
712: 
714:   *flg = SAME_NONZERO_PATTERN;
715:   MatCopy(user->M,*J,*flg);
716:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
717:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
718:   return(0);
719: }
722: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
723: {
725:   PetscScalar    ***l,***u;
726:   PetscInt       xs,xm,ys,ym;
727:   PetscInt       j,i;
728: 
730:   DMDAVecGetArrayDOF(da,xl,&l);
731:   DMDAVecGetArrayDOF(da,xu,&u);
732: 
733:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
734: 
735:   for (j=ys; j<ys+ym; j++) {
736:     for(i=xs; i < xs+xm;i++) {
737:       l[j][i][0] = -SNES_VI_INF;
738:       l[j][i][1] = 0.0;
739:       l[j][i][2] = -SNES_VI_INF;
740:       l[j][i][3] = 0.0;
741:       l[j][i][4] = 0.0;
742:       u[j][i][0] = SNES_VI_INF;
743:       u[j][i][1] = 1.0;
744:       u[j][i][2] = SNES_VI_INF;
745:       u[j][i][3] = 1.0;
746:       u[j][i][4] = 1.0;
747:     }
748:   }
749: 
750: 
751:   DMDAVecRestoreArrayDOF(da,xl,&l);
752:   DMDAVecRestoreArrayDOF(da,xu,&u);
753:   return(0);
754: }

758: PetscErrorCode GetParams(AppCtx* user)
759: {
761:   PetscBool      flg,graphicsfile = PETSC_FALSE;
762: 
764: 
765:   /* Set default parameters */
766:   user->xmin = 0.0; user->xmax = 128.0;
767:   user->ymin = 0.0; user->ymax = 128.0;
768:   user->Dv    = 1.0;
769:   user->Di    = 4.0;
770:   user->Evf   = 0.8;
771:   user->Eif   = 1.2;
772:   user->A     = 1.0;
773:   user->kBT   = 0.11;
774:   user->kav   = 1.0;
775:   user->kai   = 1.0;
776:   user->kaeta = 1.0;
777:   user->Rsurf = 10.0;
778:   user->Rbulk = 1.0;
779:   user->VG    = 100.0;
780:   user->L     = 10.0;

782:   user->T          = 1.0e-2;
783:   user->dt         = 1.0e-4;
784:   user->initv      = .00069;
785:   user->initeta    = 0.0;
786:   user->degenerate = PETSC_FALSE;
787:   user->maxevents  = 1000;
788:   user->graphics   = PETSC_TRUE;

790:   /* twodomain modeling */
791:   user->twodomain = PETSC_FALSE;
792:   user->Svr       = 0.5;
793:   user->Sir       = 0.5;
794:   user->cv_eq     = 6.9e-4;
795:   user->ci_eq     = 6.9e-4;
796:   /* void growth */
797:   user->voidgrowth = PETSC_FALSE;

799:   user->radiation = PETSC_FALSE;

801:   /* degenerate mobility */
802:   user->smallnumber = 1.0e-3;
803:   PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Coupled Cahn-Hillard/Allen-Cahn Equations","Phasefield");
804:     PetscOptionsReal("-Dv","Vacancy Diffusivity\n","None",user->Dv,&user->Dv,&flg);
805:     PetscOptionsReal("-Di","Interstitial Diffusivity\n","None",user->Di,&user->Di,&flg);
806:     PetscOptionsReal("-Evf","Vacancy Formation Energy\n","None",user->Evf,&user->Evf,&flg);
807:     PetscOptionsReal("-Eif","Interstitial Formation energy\n","None",user->Eif,&user->Eif,&flg);
808:     PetscOptionsReal("-A","???","None",user->A,&user->A,&flg);
809:     PetscOptionsReal("-kBT","Boltzmann's Constant times the Absolute Temperature","None",user->kBT,&user->kBT,&flg);
810:     PetscOptionsReal("-kav","???","None",user->kav,&user->kav,&flg);
811:     PetscOptionsReal("-kai","???","None",user->kai,&user->kai,&flg);
812:     PetscOptionsReal("-kaeta","???","None",user->kaeta,&user->kaeta,&flg);
813:     PetscOptionsReal("-Rsurf","???","None",user->Rsurf,&user->Rsurf,&flg);
814:     PetscOptionsReal("-Rbulk","???","None",user->Rbulk,&user->Rbulk,&flg);
815:     PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
816:     PetscOptionsReal("-L","???","None",user->L,&user->L,&flg);

818:     PetscOptionsReal("-initv","Initial solution of Cv and Ci","None",user->initv,&user->initv,&flg);
819:     PetscOptionsReal("-initeta","Initial solution of Eta","None",user->initeta,&user->initeta,&flg);
820:     PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
821:     PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
822:     PetscOptionsBool("-twodomain","Run two domain model\n","None",user->twodomain,&user->twodomain,&flg);
823:     PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
824:     PetscOptionsBool("-radiation","Use initial conditions for void growth\n","None",user->radiation,&user->radiation,&flg);
825:     PetscOptionsReal("-xmin","Lower X coordinate of domain\n","None",user->xmin,&user->xmin,&flg);
826:     PetscOptionsReal("-xmax","Upper X coordinate of domain\n","None",user->xmax,&user->xmax,&flg);
827:     PetscOptionsReal("-T","Total runtime\n","None",user->T,&user->T,&flg);
828:     PetscOptionsReal("-dt","Time step\n","None",user->dt,&user->dt,&flg);
829:     user->dtevent = user->dt;
830:     PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);
831:     PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);

833:     PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
834:     PetscOptionsBool("-graphicsfile","Save solution at each timestep\n","None",graphicsfile,&graphicsfile,&flg);
835:     if (graphicsfile) {
836:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ex61.data",FILE_MODE_WRITE,&user->graphicsfile);
837:     }
838:   PetscOptionsEnd();
839:   return(0);
840:  }


845: PetscErrorCode SetUpMatrices(AppCtx* user)
846: {
847:   PetscErrorCode    ierr;
848:   PetscInt          nele,nen,i,n;
849:   const PetscInt    *ele;
850:   PetscScalar       dt=user->dt,hx,hy;
851: 
852:   PetscInt          idx[3];
853:   PetscScalar       eM_0[3][3],eM_2_even[3][3],eM_2_odd[3][3];
854:   PetscScalar       cv_sum, ci_sum;
855:   Mat               M=user->M;
856:   Mat               M_0=user->M_0;
857:   PetscInt          Mda=user->Mda, Nda=user->Nda;
858:   PetscScalar       *cv_p,*ci_p;
859:   /* newly added */
860:   Vec               cvlocal,cilocal;

863: 
864:   /*  MatSetOption(M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
865:    MatSetOption(M_0,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);*/

867:   /* new stuff */
868:   DMGetLocalVector(user->da2,&cvlocal);
869:   DMGetLocalVector(user->da2,&cilocal);
870:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
871:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
872:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
873:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
874:   /* old stuff */
875:   /*
876:   VecGetArray(user->cv,&cv_p);
877:   VecGetArray(user->ci,&ci_p);
878:    */
879:   /* new stuff */
880:   VecGetArray(cvlocal,&cv_p);
881:   VecGetArray(cilocal,&ci_p);

883:   MatGetLocalSize(M,&n,PETSC_NULL);
884:   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);
885:   hx = (user->xmax-user->xmin)/Mda;
886:   hy = (user->ymax-user->ymin)/Nda;

888:   eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hy/12.0;
889:   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;
890: 
891:   eM_2_odd[0][0] = 1.0;
892:   eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
893:   eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
894:   eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;

896:   eM_2_even[0][0] = 1.0;
897:   eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
898:   eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
899:   eM_2_even[1][2] = eM_2_even[2][1] = 0.0;

901:   /*  eM_2_even[1][1] = 1.0;
902:   eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
903:   eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
904:   eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
905:    */

907:   //  for(k=0;k < Mda*Nda*2;k++) {
908:   DMDAGetElements(user->da1,&nele,&nen,&ele);
909:   for (i=0; i < nele; i++) {
910:     /*
911:     idx[0] = connect[k*3];
912:     idx[1] = connect[k*3+1];
913:     idx[2] = connect[k*3+2];
914:      */
915:     idx[0] = ele[3*i];
916:     idx[1] = ele[3*i+1];
917:     idx[2] = ele[3*i+2];

919:     PetscInt    row,cols[6],r,row_M_0,cols3[3];
920:     PetscScalar vals[6],vals_M_0[3],vals3[3];
921: 
922:     for(r=0;r<3;r++) {
923:       //row_M_0 = connect[k*3+r];
924:       row_M_0 = idx[r];

926:       vals_M_0[0]=eM_0[r][0];
927:       vals_M_0[1]=eM_0[r][1];
928:       vals_M_0[2]=eM_0[r][2];
929: 
930: 
931:       MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
932: 
933:       if (user->degenerate) {
934:         cv_sum = (cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
935:         ci_sum = (ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
936:       } else {
937:         cv_sum = user->initv*user->Dv/user->kBT;
938:         ci_sum = user->initv*user->Di/user->kBT;
939:       }
940: 

942: 
943:         row = 5*idx[r];
944:         cols[0] = 5*idx[0];     vals[0] = dt*eM_2_odd[r][0]*cv_sum;
945:         cols[1] = 5*idx[1];     vals[1] = dt*eM_2_odd[r][1]*cv_sum;
946:         cols[2] = 5*idx[2];     vals[2] = dt*eM_2_odd[r][2]*cv_sum;
947:         cols[3] = 5*idx[0]+1;   vals[3] = eM_0[r][0];
948:         cols[4] = 5*idx[1]+1;   vals[4] = eM_0[r][1];
949:         cols[5] = 5*idx[2]+1;   vals[5] = eM_0[r][2];

951: 
952:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
953: 
954: 
955:         row = 5*idx[r]+1;
956:         cols[0] = 5*idx[0];     vals[0] = -1.0*eM_0[r][0];
957:         cols[1] = 5*idx[1];     vals[1] = -1.0*eM_0[r][1];
958:         cols[2] = 5*idx[2];     vals[2] = -1.0*eM_0[r][2];
959:         cols[3] = 5*idx[0]+1;   vals[3] =  user->kav*eM_2_odd[r][0];
960:         cols[4] = 5*idx[1]+1;   vals[4] =  user->kav*eM_2_odd[r][1];
961:         cols[5] = 5*idx[2]+1;   vals[5] =  user->kav*eM_2_odd[r][2];
962: 
963:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
964: 
965:         row = 5*idx[r]+2;
966:         cols[0] = 5*idx[0]+2;   vals[0] =  dt*eM_2_odd[r][0]*ci_sum;
967:         cols[1] = 5*idx[1]+2;   vals[1] =  dt*eM_2_odd[r][1]*ci_sum;
968:         cols[2] = 5*idx[2]+2;   vals[2] =  dt*eM_2_odd[r][2]*ci_sum;
969:         cols[3] = 5*idx[0]+3;   vals[3] =  eM_0[r][0];
970:         cols[4] = 5*idx[1]+3;   vals[4] =  eM_0[r][1];
971:         cols[5] = 5*idx[2]+3;   vals[5] =  eM_0[r][2];
972: 
973:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
974: 
975: 
976:         row = 5*idx[r]+3;
977:         cols[0] = 5*idx[0]+2;   vals[0] = -1.0*eM_0[r][0];
978:         cols[1] = 5*idx[1]+2;   vals[1] = -1.0*eM_0[r][1];
979:         cols[2] = 5*idx[2]+2;   vals[2] = -1.0*eM_0[r][2];
980:         cols[3] = 5*idx[0]+3;   vals[3] =  user->kai*eM_2_odd[r][0];
981:         cols[4] = 5*idx[1]+3;   vals[4] =  user->kai*eM_2_odd[r][1];
982:         cols[5] = 5*idx[2]+3;   vals[5] =  user->kai*eM_2_odd[r][2];
983: 
984:         MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
985: 
986: 
987:         row = 5*idx[r]+4;
988:         /*
989:         cols3[0] = 5*idx[0]+4;   vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_odd[r][0];
990:         cols3[1] = 5*idx[1]+4;   vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_odd[r][1];
991:         cols3[2] = 5*idx[2]+4;   vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_odd[r][2];
992:          */
993:         cols3[0] = 5*idx[0]+4;   vals3[0] = (eM_0[r][0]/dt + user->L*user->kaeta*eM_2_odd[r][0]);
994:         cols3[1] = 5*idx[1]+4;   vals3[1] = (eM_0[r][1]/dt + user->L*user->kaeta*eM_2_odd[r][1]);
995:         cols3[2] = 5*idx[2]+4;   vals3[2] = (eM_0[r][2]/dt + user->L*user->kaeta*eM_2_odd[r][2]);
996: 
997:         MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);

999:     }
1000:   }

1002:   /* new */
1003:   VecRestoreArray(cvlocal,&cv_p);
1004:   VecRestoreArray(cilocal,&ci_p);
1005:   DMRestoreLocalVector(user->da2,&cvlocal);
1006:   DMRestoreLocalVector(user->da2,&cilocal);
1007:   /* old */
1008:   /*
1009:   VecRestoreArray(user->cv,&cv_p);
1010:   VecRestoreArray(user->ci,&ci_p);
1011:    */
1012:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
1013:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
1014: 
1015:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1016:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1017: 
1018:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
1019: 

1021:   return(0);
1022: }


1027: PetscErrorCode UpdateMatrices(AppCtx* user)
1028: {
1029:   PetscErrorCode    ierr;
1030:   PetscInt          i,n,Mda,Nda,nele,nen;
1031:   const PetscInt    *ele;
1032: 
1033:   PetscInt          idx[3];
1034:   PetscScalar       eM_2_odd[3][3],eM_2_even[3][3],h,dt=user->dt;
1035:   Mat               M=user->M;
1036:   PetscScalar       *cv_p,*ci_p,cv_sum,ci_sum;
1037:   /* newly added */
1038:   Vec               cvlocal,cilocal;

1041: 
1042: 
1043:   MatGetLocalSize(M,&n,PETSC_NULL);
1044: 
1045:   /* new stuff */
1046:   DMGetLocalVector(user->da2,&cvlocal);
1047:   DMGetLocalVector(user->da2,&cilocal);
1048:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
1049:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
1050:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
1051:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
1052:   /* new stuff */
1053:   VecGetArray(cvlocal,&cv_p);
1054:   VecGetArray(cilocal,&ci_p);
1055: 
1056:   /* old stuff */
1057:   /*
1058:   VecGetArray(user->cv,&cv_p);
1059:   VecGetArray(user->ci,&ci_p);
1060:    */
1061:   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);

1063: 
1064: 
1065:   h = (user->xmax-user->xmin)/Mda;

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

1069:   for(i=0; i < nele; i++) {
1070:     /*
1071:     idx[0] = connect[k*3];
1072:     idx[1] = connect[k*3+1];
1073:     idx[2] = connect[k*3+2];
1074:      */
1075:     idx[0] = ele[3*i];
1076:     idx[1] = ele[3*i+1];
1077:     idx[2] = ele[3*i+2];

1079:     PetscInt r,row,cols[3];
1080:     PetscScalar vals[3];
1081:     for (r=0;r<3;r++) {
1082:       row = 5*idx[r];
1083:       cols[0] = 5*idx[0];     vals[0] = 0.0;
1084:       cols[1] = 5*idx[1];     vals[1] = 0.0;
1085:       cols[2] = 5*idx[2];     vals[2] = 0.0;

1087:       /* Insert values in matrix M for 1st dof */
1088:       MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1089: 
1090:       row = 5*idx[r]+2;
1091:       cols[0] = 5*idx[0]+2;   vals[0] = 0.0;
1092:       cols[1] = 5*idx[1]+2;   vals[1] = 0.0;
1093:       cols[2] = 5*idx[2]+2;   vals[2] = 0.0;

1095:       /* Insert values in matrix M for 3nd dof */
1096:       MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1097:     }
1098:   }
1099: 
1100:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1101:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

1103:   eM_2_odd[0][0] = 1.0;
1104:   eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
1105:   eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
1106:   eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;

1108:   eM_2_even[0][0] = 1.0;
1109:   eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
1110:   eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
1111:   eM_2_even[1][2] = eM_2_even[2][1] = 0.0;

1113:   /*
1114:   eM_2_even[1][1] = 1.0;
1115:   eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
1116:   eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
1117:   eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
1118:    */
1119: 
1120:   /* Get local element info */
1121:   //for(k=0;k < Mda*Nda*2;k++) {
1122:   for (i=0; i < nele; i++) {
1123:     /*
1124:       idx[0] = connect[k*3];
1125:       idx[1] = connect[k*3+1];
1126:       idx[2] = connect[k*3+2];
1127:      */
1128:     idx[0] = ele[3*i];
1129:     idx[1] = ele[3*i+1];
1130:     idx[2] = ele[3*i+2];

1132:       PetscInt    row,cols[3],r;
1133:       PetscScalar vals[3];
1134: 
1135:       for(r=0;r<3;r++) {
1136: 
1137:       if (user->degenerate) {
1138:         printf("smallnumber = %14.12f\n",user->smallnumber);
1139:         cv_sum = (user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
1140:         ci_sum = (user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
1141:       } else {
1142:         cv_sum = user->initv*user->Dv/(user->kBT);
1143:         ci_sum = user->initv*user->Di/user->kBT;
1144:       }

1146: 
1147: 
1148:                 row = 5*idx[r];
1149:                 cols[0] = 5*idx[0];     vals[0] = dt*eM_2_odd[r][0]*cv_sum;
1150:                 cols[1] = 5*idx[1];     vals[1] = dt*eM_2_odd[r][1]*cv_sum;
1151:                 cols[2] = 5*idx[2];     vals[2] = dt*eM_2_odd[r][2]*cv_sum;
1152: 
1153:                 /* Insert values in matrix M for 1st dof */
1154:                 MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
1155: 
1156: 
1157:                 row = 5*idx[r]+2;
1158:                 cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2_odd[r][0]*ci_sum;
1159:                 cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2_odd[r][1]*ci_sum;
1160:                 cols[2] = 5*idx[2]+2;   vals[2] = dt*eM_2_odd[r][2]*ci_sum;
1161: 
1162:                 MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);

1164: 
1165: 
1166:         }
1167: 
1168:     }

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

1181:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1182:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);


1185:   return(0);
1186: }



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


1204:   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);
1205:   DMDAGetGhostedCoordinates(user->da2,&coords);
1206:   VecGetArrayRead(coords,&_coords);

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


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

1219:     x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
1220:     x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
1221:     x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

1223:     //printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);
1224:     //printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);
1225: 
1226:     PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
1227:     PetscInt    k;

1229:     xc1 = user->xmin;
1230:     xc2 = xmid;

1232:     //VecSet(user->phi1,0.0);
1233:     for (k=0;k < 3; k++) {
1234:       if (x[k]-xqu > 0) {
1235:         s1 = (x[k] - xqu);
1236:       } else {
1237:         s1 = -(x[k] - xqu);
1238:       }
1239:       if (x[k] - xc1 > 0) {
1240:         dist1 = (x[k] - xc1);
1241:       } else {
1242:         dist1 = -(x[k] - xc1);
1243:       }
1244:       if (x[k] - xc2 > 0) {
1245:         dist2 = (x[k] - xc2);
1246:       } else {
1247:         dist2 = -(x[k] - xc2);
1248:       }
1249:       if (dist1 <= 0.5*lambda) {
1250:         r = (x[k]-xc1)/(0.5*lambda);
1251:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1252:         vals1[k] = hhr;
1253:       }
1254:       else if (dist2 <= 0.5*lambda) {
1255:         r = (x[k]-xc2)/(0.5*lambda);
1256:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1257:         vals1[k] = 1.0 - hhr;
1258:       }
1259:       else if (s1 <= xqu - 2.0*h) {
1260:         vals1[k] = 1.0;
1261:       }
1262: 
1263:       //else if ( abs(x[k]-(user->xmax-h)) < 0.1*h ) {
1264:       else if ( (user->xmax-h)-x[k] < 0.1*h ) {
1265:         vals1[k] = .15625;
1266:        }
1267:       else {
1268:         vals1[k] = 0.0;
1269:       }
1270:     }
1271: 
1272:     VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);

1274:     xc1 = xmid;
1275:     xc2 = user->xmax;

1277:     //VecSet(user->phi2,0.0);
1278:     for (k=0;k < 3; k++) {
1279:       /*
1280:       s1 = abs(x[k] - (xqu+xmid));
1281:       dist1 = abs(x[k] - xc1);
1282:       dist2 = abs(x[k] - xc2);
1283:        */
1284:       if (x[k]-(xqu + xmid) > 0) {
1285:         s1 = (x[k] - (xqu + xmid));
1286:       } else {
1287:         s1 = -(x[k] - (xqu + xmid));
1288:       }
1289:       if (x[k] - xc1 > 0) {
1290:         dist1 = (x[k] - xc1);
1291:       } else {
1292:         dist1 = -(x[k] - xc1);
1293:       }
1294:       if (x[k] - xc2 > 0) {
1295:         dist2 = (x[k] - xc2);
1296:       } else {
1297:         dist2 = -(x[k] - xc2);
1298:       }
1299: 
1300:       if (dist1 <= 0.5*lambda) {
1301:         r = (x[k]-xc1)/(0.5*lambda);
1302:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1303:         vals2[k] = hhr;
1304:       }
1305:       else if (dist2 <= 0.5*lambda) {
1306:         r = -(x[k]-xc2)/(0.5*lambda);
1307:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1308:         vals2[k] = hhr;
1309:       }
1310:       else if (s1 <= xqu - 2.0*h) {
1311:         vals2[k] = 1.0;
1312:       }
1313: 
1314:       else if ( x[k]-(user->xmin) < 0.1*h ) {
1315:         vals2[k] = 0.5;
1316:       }
1317: 
1318: 
1319:       else if ( (x[k]-(user->xmin+h)) < 0.1*h ) {
1320:         vals2[k] = .15625;
1321:       }
1322: 
1323:       else {
1324:         vals2[k] = 0.0;
1325:       }
1326: 
1327:     }

1329:     VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
1330:     /*
1331:     for (k=0;k < 3; k++) {
1332:       vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
1333:     }
1334:      */
1335:     //VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);
1336: 
1337:   }
1338: 
1339:   VecAssemblyBegin(user->phi1);
1340:   VecAssemblyEnd(user->phi1);
1341:   VecAssemblyBegin(user->phi2);
1342:   VecAssemblyEnd(user->phi2);
1343:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);

1345:   VecView(user->phi1,view);
1346:   VecView(user->phi2,view);

1348: 
1349:   //VecView(user->phi1,0);
1350:   //VecView(user->phi2,0);
1351: 
1352:   VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1353:   VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1354:   VecView(user->phi1,view);
1355:   VecView(user->phi2,view);

1357:   VecCopy(user->phi1,user->Phi2D_V);
1358:   VecAXPY(user->Phi2D_V,1.0,user->phi2);
1359:   //VecView(user->Phi2D_V,0);

1361:   VecView(user->Phi2D_V,view);
1362:   PetscViewerDestroy(&view);
1363:   //  VecNorm(user->Phi2D_V,NORM_INFINITY,&max1);
1364:   //VecMin(user->Phi2D_V,&loc1,&min1);
1365:   //printf("norm phi = %f, min phi = %f\n",max1,min1);

1367:   return(0);
1368: 
1369: }