Commits

Nathan Collier committed 995248a

working changes to PFC

Comments (0)

Files changed (1)

demo/PhaseFieldCrystal2D.c

 
 typedef struct { 
   IGA iga;
-  PetscReal theta,cbar,alpha;
-  PetscReal L0,lambda,Sprev[3];
+  PetscReal cbar,D,k,Eps,g;
+  PetscReal L0,Sprev[3];
 } AppCtx;
 
 #undef __FUNCT__
 #define __FUNCT__ "GinzburgLandauFreeEnergy"
 PetscScalar GinzburgLandauFreeEnergy(PetscReal c,PetscReal cx,PetscReal cy,AppCtx *user)
 {
-  PetscReal E,omc = 1.0-c;
-  E=c*log(c)+omc*log(omc)+2.0*user->theta*c*omc+user->theta/3.0/user->alpha*fabs(cx*cx+cy*cy);
+  PetscReal E;
+  E=-user->Eps/2.0*c*c-user->g/3.0*c*c*c+c*c*c/4.0+user->D/2.0*((cx*cx+cy*cy)*(cx*cx+cy*cy)-2.0*user->k*user->k*fabs(cx*cx+cy*cy)+user->k*user->k*user->k*user->k*c*c);
   return E;
 }
 
   PetscFunctionReturn(0);
 }
 
-#undef  __FUNCT__
-#define __FUNCT__ "Mobility"
-void Mobility(AppCtx *user,PetscReal c,PetscReal *M,PetscReal *dM,PetscReal *d2M)
-{
-  if (M)   *M   = c*(1-c);
-  if (dM)  *dM  = 1-2*c;
-  if (d2M) *d2M = -2;
-}
+//#undef  __FUNCT__
+//#define __FUNCT__ "Mobility"
+//void Mobility(AppCtx *user,PetscReal c,PetscReal *M,PetscReal *dM,PetscReal *d2M)
+//{
+//  if (M)   *M   = c*(1-c);
+//  if (dM)  *dM  = 1-2*c;
+//  if (d2M) *d2M = -2;
+//}
 
 #undef __FUNCT__
 #define __FUNCT__ "StatsMonitor"
   PetscFunctionReturn(0);
 }
 
-#undef  __FUNCT__
-#define __FUNCT__ "ChemicalPotential"
-void ChemicalPotential(AppCtx *user,PetscReal c,PetscReal *mu,PetscReal *dmu,PetscReal *d2mu)
-{
-  if (mu) {
-    (*mu)  = 0.5/user->theta*log(c/(1-c))+1-2*c;
-    (*mu) *= user->L0*user->L0/user->lambda;
-  }
-  if (dmu) {
-    (*dmu)  = 0.5/user->theta*1.0/(c*(1-c)) - 2;
-    (*dmu) *= user->L0*user->L0/user->lambda;
-  }
-  if (d2mu) {
-    (*d2mu)  = -0.5/user->theta*(1-2*c)/(c*c*(1-c)*(1-c));
-    (*d2mu) *= user->L0*user->L0/user->lambda;
-  }
-}
+//#undef  __FUNCT__
+//#define __FUNCT__ "ChemicalPotential"
+//void ChemicalPotential(AppCtx *user,PetscReal c,PetscReal *mu,PetscReal *dmu,PetscReal *d2mu)
+//{
+//  if (mu) {
+//   (*mu)  = 0.5/user->theta*log(c/(1-c))+1-2*c;
+//    (*mu) *= user->L0*user->L0/user->lambda;
+//  }
+//  if (dmu) {
+//    (*dmu)  = 0.5/user->theta*1.0/(c*(1-c)) - 2;
+//    (*dmu) *= user->L0*user->L0/user->lambda;
+//  }
+//  if (d2mu) {
+//    (*d2mu)  = -0.5/user->theta*(1-2*c)/(c*c*(1-c)*(1-c));
+//    (*d2mu) *= user->L0*user->L0/user->lambda;
+//  }
+//}
 
 #undef  __FUNCT__
 #define __FUNCT__ "Residual"
   IGAPointGetValue(p,V,&c_t);
   IGAPointGetValue(p,U,&c);
 
-  PetscReal M,dM;
-  Mobility(user,c,&M,&dM,NULL);
-  PetscReal dmu;
-  ChemicalPotential(user,c,NULL,&dmu,NULL);
+  //PetscReal M,dM;
+  //Mobility(user,c,&M,&dM,NULL);
+  //PetscReal dmu;
+  //ChemicalPotential(user,c,NULL,&dmu,NULL);
 
-  PetscScalar c1[2],c2[2][2];
+  PetscScalar c1[2],c2[2][2],c3[2][2][2];
   IGAPointGetGrad(p,U,&c1[0]);
   IGAPointGetHess(p,U,&c2[0][0]);
-  PetscScalar c_x  = c1[0],    c_y  = c1[1];
-  PetscScalar c_xx = c2[0][0], c_yy = c2[1][1];
+  IGAPointGet3rdMixedPartials(p,U,&c3[0][0][0]);
+  PetscScalar c_x   = c1[0],       c_y   = c1[1];
+  PetscScalar c_xx  = c2[0][0],    c_yy  = c2[1][1];
+  PetscScalar c_xxx = c3[0][0][0], c_yyy = c3[1][1][1];
+  PetscScalar c_xxy = c3[0][0][1], c_xyy = c3[0][1][1];
 
-  const PetscReal *N0,(*N1)[2],(*N2)[2][2];
+  const PetscReal *N0,(*N1)[2],(*N2)[2][2],(*N3)[2][2][2];
   IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
   IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+  IGAPointGetShapeFuns(p,3,(const PetscReal**)&N3);
 
   PetscInt a;
   for (a=0; a<nen; a++) {
-    PetscReal Na    = N0[a];
-    PetscReal Na_x  = N1[a][0];
-    PetscReal Na_y  = N1[a][1];
-    PetscReal Na_xx = N2[a][0][0];
-    PetscReal Na_yy = N2[a][1][1];
+    PetscReal Na     = N0[a];
+    PetscReal Na_x   = N1[a][0];
+    PetscReal Na_y   = N1[a][1];
+    PetscReal Na_xx  = N2[a][0][0];
+    PetscReal Na_yy  = N2[a][1][1];
+    PetscReal Na_xxx = N3[a][0][0][0];
+    PetscReal Na_yyy = N3[a][1][1][1];
+    PetscReal Na_xxy = N3[a][0][0][1];
+    PetscReal Na_xyy = N3[a][0][1][1];
+
     /* ----- */
     PetscScalar Ra  = 0;
     // Na * c_t
     Ra += Na * c_t; 
-    // grad(Na) . ((M*dmu + dM*del2(c))) grad(C)
-    PetscScalar t1 = M*dmu + dM*(c_xx+c_yy);
+    // grad(Na) . ((-\epsilon\phi-g\epsilon^2+\phi^3)) grad(c)
+    PetscScalar t1 = -user->Eps*c - user->g*c*c +c*c*c;
     Ra += Na_x * t1 * c_x;
     Ra += Na_y * t1 * c_y;
-    // del2(Na) * M * del2(c)
-    Ra += (Na_xx+Na_yy) * M * (c_xx+c_yy);
+    // -2Dk^2 . del2(Na)  . del2(c)
+    Ra += -2.0*user->D*user->k*user->k * (Na_xx+Na_yy) * (c_xx+c_yy);
+    // D . del3(Na) . del3(c)
+    Ra += user->D * (Na_yyy * (c_yyy + c_xxy) + Na_xxx * (c_xxx + c_xyy) + Na_xyy * (c_xxx + c_xyy) + Na_xxy * (c_yyy + c_xxy));
     /* ----- */
     R[a] = Ra;
   }
   IGAPointGetValue(p,V,&c_t);
   IGAPointGetValue(p,U,&c);
 
-  PetscReal M,dM,d2M;
-  Mobility(user,c,&M,&dM,&d2M);
-  PetscReal dmu,d2mu;
-  ChemicalPotential(user,c,NULL,&dmu,&d2mu);
+  //PetscReal M,dM,d2M;
+  //Mobility(user,c,&M,&dM,&d2M);
+  // PetscReal dmu,d2mu;
+  // ChemicalPotential(user,c,NULL,&dmu,&d2mu);
 
-  PetscScalar c1[2],c2[2][2];
+  PetscScalar c1[2],c2[2][2],c3[2][2][2];
   IGAPointGetGrad(p,U,&c1[0]);
   IGAPointGetHess(p,U,&c2[0][0]);
-  PetscScalar c_x  = c1[0],    c_y  = c1[1];
-  PetscScalar c_xx = c2[0][0], c_yy = c2[1][1];
+  IGAPointGet3rdMixedPartials(p,U,&c3[0][0][0]);
+  PetscScalar c_x   = c1[0],       c_y   = c1[1];
+  PetscScalar c_xx  = c2[0][0],    c_yy  = c2[1][1];
+  PetscScalar c_xxx = c3[0][0][0], c_yyy = c3[1][1][1];
+  PetscScalar c_xxy = c3[0][0][1], c_xyy = c3[0][1][1];
 
-  const PetscReal *N0,(*N1)[2],(*N2)[2][2];
+  const PetscReal *N0,(*N1)[2],(*N2)[2][2],(*N3)[2][2][2];
   IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
   IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+  IGAPointGetShapeFuns(p,3,(const PetscReal**)&N3);
 
   PetscInt a,b;
   for (a=0; a<nen; a++) {
-    PetscReal Na    = N0[a];
-    PetscReal Na_x  = N1[a][0];
-    PetscReal Na_y  = N1[a][1];
-    PetscReal Na_xx = N2[a][0][0];
-    PetscReal Na_yy = N2[a][1][1];
+    PetscReal Na     = N0[a];
+    PetscReal Na_x   = N1[a][0];
+    PetscReal Na_y   = N1[a][1];
+    PetscReal Na_xx  = N2[a][0][0];
+    PetscReal Na_yy  = N2[a][1][1];
+    PetscReal Na_xxx = N3[a][0][0][0];
+    PetscReal Na_yyy = N3[a][1][1][1];
+    PetscReal Na_xxy = N3[a][0][0][1];
+    PetscReal Na_xyy = N3[a][0][1][1];
+
     for (b=0; b<nen; b++) {
-      PetscReal Nb    = N0[b];
-      PetscReal Nb_x  = N1[b][0];
-      PetscReal Nb_y  = N1[b][1];
-      PetscReal Nb_xx = N2[b][0][0];
-      PetscReal Nb_yy = N2[b][1][1];
+      PetscReal Nb     = N0[b];
+      PetscReal Nb_x   = N1[b][0];
+      PetscReal Nb_y   = N1[b][1];
+      PetscReal Nb_xx  = N2[b][0][0];
+      PetscReal Nb_yy  = N2[b][1][1];
+      PetscReal Nb_xxx = N3[b][0][0][0];
+      PetscReal Nb_yyy = N3[b][1][1][1];
+      PetscReal Nb_xxy = N3[b][0][0][1];
+      PetscReal Nb_xyy = N3[b][0][1][1];
+
       /* ----- */
       PetscScalar Kab = 0;
       // shift*Na*Nb
       Kab += shift*Na*Nb;
-      // grad(Na) . (M*dmu+dM*del2(c)) grad(Nb)
-      PetscScalar t1 = M*dmu + dM*(c_xx+c_yy);
+      // grad(Na) . (-Eps*c-g*c^2+c^3) grad(Nb)
+      PetscScalar t1 = -user->Eps*c - user->g*c*c +c*c*c;
       Kab += Na_x * t1 * Nb_x;
       Kab += Na_y * t1 * Nb_y;
-      // grad(Na) . ((dM*dmu+M*d2mu+d2M*del2(c))*Nb + dM*del2(Nb)) grad(C)
-      PetscScalar t2 = (dM*dmu+M*d2mu+d2M*(c_xx+c_yy))*Nb + dM*(Nb_xx+Nb_yy);
+      // grad(Na) . (deriv(t1)) grad(c)
+      PetscScalar t2 = (-user->Eps*Nb-user->g*Nb*Nb+Nb*Nb*Nb);
       Kab += Na_x * t2 * c_x;
       Kab += Na_y * t2 * c_y;
-      // del2(Na) * ((dM*del2(c)*Nb + M*del2(Nb))
-      Kab += (Na_xx+Na_yy) * (dM*(c_xx+c_yy)*Nb + M*(Nb_xx+Nb_yy));
+      // -2Dk^2 * del2(Na) * del2(Nb)
+      Kab += -2.0*user->D*user->k*user->k*(Na_xx+Na_yy) * (Nb_xx+Nb_yy);
+      // D * del3(Na) * del3(Nb)
+      Kab += user->D*(Na_yyy * (Nb_yyy + Nb_xxy) + Na_xxx * (Nb_xxx + Nb_xyy) + Na_xyy * (Nb_xxx + Nb_xyy) + Na_xxy * (Nb_yyy + Nb_xxy));
       /* ----- */
       K[a*nen+b] = Kab;
     }
     PetscRandom rctx;    
     ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
     ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
-    ierr = PetscRandomSetInterval(rctx,user->cbar-0.05,user->cbar+0.05);CHKERRQ(ierr); 
+    ierr = PetscRandomSetInterval(rctx,user->cbar-0.000001,user->cbar+0.0000001);CHKERRQ(ierr); 
     ierr = PetscRandomSeed(rctx);CHKERRQ(ierr);
     ierr = VecSetRandom(C,rctx);CHKERRQ(ierr); 
     ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 
 
   /* Define simulation specific parameters */
   AppCtx user;
-  user.cbar  = 0.63;   /* average concentration */
-  user.alpha = 3000.0; /* thickess interface parameter */
-  user.theta = 1.5;    /* temperature/critical temperature */
-  user.L0    = 1.0;    /* length scale */
+  user.cbar  = 0.285;   /* average concentration */
+  user.D     = 1.0;     /* thickess interface parameter */
+  user.k     = 1.0;     /* thickess interface parameter */
+  user.Eps   = 0.25;    /* thickess interface parameter */
+  user.g     = 0.0;     /* thickess interface parameter */
+  user.L0    = 800.0;   /* length scale */
   user.Sprev[0] = user.Sprev[1] = user.Sprev[2] = 1.0e20; 
 
   /* Set discretization options */
   ierr = PetscOptionsString("-ch_initial","Load initial solution from file",__FILE__,initial,initial,sizeof(initial),PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-ch_output","Enable output files",__FILE__,output,&output,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-ch_monitor","Compute and show statistics of solution",__FILE__,monitor,&monitor,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-ch_cbar","Initial average concentration",__FILE__,user.cbar,&user.cbar,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-ch_alpha","Characteristic parameter",__FILE__,user.alpha,&user.alpha,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-ch_cbar","Initial atomistic density field",__FILE__,user.cbar,&user.cbar,PETSC_NULL);CHKERRQ(ierr);
+  //ierr = PetscOptionsReal("-ch_alpha","Characteristic parameter",__FILE__,user.alpha,&user.alpha,PETSC_NULL);CHKERRQ(ierr);
+ierr = PetscOptionsReal("-ch_g","Physical parameter",__FILE__,user.g,&user.g,PETSC_NULL);CHKERRQ(ierr);
+ierr = PetscOptionsReal("-ch_k","Positive number",__FILE__,user.k,&user.k,PETSC_NULL);CHKERRQ(ierr);
+ierr = PetscOptionsReal("-ch_Eps","Physical parameter",__FILE__,user.Eps,&user.Eps,PETSC_NULL);CHKERRQ(ierr);
+ierr = PetscOptionsReal("-ch_D","Positive number",__FILE__,user.D,&user.D,PETSC_NULL);CHKERRQ(ierr);
+
+
+
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
   if (C == PETSC_DECIDE) C = p-1;
 
-  user.lambda = 1.0/N/N; /* mesh size parameter */
+  // user.lambda = 1.0/N/N; /* mesh size parameter */
   
   if (p < 2 || C < 1) /* Problem requires a p>=2 C1 basis */
     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,
   ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
   ierr = IGAAxisSetPeriodic(axis0,PETSC_TRUE);CHKERRQ(ierr);
   ierr = IGAAxisSetDegree(axis0,p);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,N,0.0,1.0,C);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,N,0.0,user.L0,C);CHKERRQ(ierr);
   IGAAxis axis1;
   ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
   ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
-  ierr = TSSetDuration(ts,10000,1.0);CHKERRQ(ierr);
-  ierr = TSSetTimeStep(ts,1e-10);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,10000000,2000.0);CHKERRQ(ierr);
+  ierr = TSSetTimeStep(ts,1e-1);CHKERRQ(ierr);
 
   ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
-  ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
-  ierr = TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,PETSC_NULL);CHKERRQ(ierr); 
+  ierr = TSAlphaSetRadius(ts,1.0);CHKERRQ(ierr);
+  //ierr = TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,PETSC_NULL);CHKERRQ(ierr); 
 
   if (monitor) {
     user.iga = iga;
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.