Commits

Vijay Mahadevan committed 93a53a1

Various fixes for the gray-scott model problem.

Comments (0)

Files changed (1)

src/dm/impls/moab/examples/tests/ex5.c

 static char help[] = "Solve the time-dependent, coupled nonlinear diffusion-reaction equations. \n\n";
 
 #include <petscts.h>
+#include <petscdt.h>
 #include <petscdmmoab.h>
 
 // MOAB includes:
 #error You must have MOAB for this example. Reconfigure using --download-moab
 #endif
 
-//#define LOCAL_ASSEMBLY
+#define LOCAL_ASSEMBLY
 #define SUBPROBLEM 1
 
 typedef struct {
   PetscReal psi;                /* neutron inverse-velocity - controls rate of diffusion   */
   PetscReal tfinal;
   PetscReal delxyz;
+  PetscReal *quadx,*quadw;
 } Data;
 
 
   TS             ts;                         /* time integrator */
   Vec            x;                        /* solution vector */
   Mat            A;
-  PetscInt       steps,maxsteps,nonlinits,linits,snesfails,rejects;  /* iterations for convergence */
+  PetscInt       steps,maxsteps,nonlinits,linits,snesfails,rejects,nquad;  /* iterations for convergence */
   PetscErrorCode ierr;
   DM             dm;
   Data           user;
   user.n      = 2;
 #if (SUBPROBLEM==1)
   // http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/
-  user.diffusion.U = 5e-4;
-  user.diffusion.V = 1e-4;
+  /*  Modeling Morphogenesis with Reaction-Diffusion Equations using Galerkin Spectral Methods
+      USNA, Trident Project
+      Benjamin M. Heineike
+      Project advisers:
+        Professor Reza Malek-Madani
+        Professor Sonia Garcia
+  */
+  user.diffusion.U = 1e-4;
+  user.diffusion.V = 1e-5;
   user.cpld_diffusion = 0.0;
-  user.eta    = 0.01;
-  user.psi    = 0.004;
+  user.eta    = 0.0154;
+  user.psi    = 0.04;
 #else
-  user.diffusion.U = 1.0e-7;
-  user.diffusion.V = 8.0e-6;
-  user.cpld_diffusion = 2.5e-6;
-  user.eta    = 0.0006;
-  user.psi    = 0.005;  
+  user.diffusion.U = 2.0e-5;
+  user.diffusion.V = 1.0e-5;
+  user.cpld_diffusion = 2.5e-6*0;
+  user.eta    = 0.0154;
+  user.psi    = 0.05;
 #endif
   maxsteps    = 10000;
   endtime     = 1.0;
+  nquad       = 2;
   use_monitor = PETSC_FALSE;
   ierr        = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the dimension-independent diffusion-reaction solver", "DMqq");
   ierr        = PetscOptionsInt("-n", "The elements in each direction", "ex5.c",user.n,&user.n, NULL);CHKERRQ(ierr);
   ierr        = PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","ex5.c",use_monitor,&use_monitor,NULL);CHKERRQ(ierr);
   ierr        = PetscOptionsInt("-tsteps","Maximum number of time steps","ex5.c",maxsteps,&maxsteps,NULL);CHKERRQ(ierr);
   ierr        = PetscOptionsReal("-endtime","Duration of temporal integration","ex5.c",endtime,&endtime,NULL);CHKERRQ(ierr);
+  ierr        = PetscOptionsInt("-nquad","Number of quadrature points in the 1-d space","ex5.c",nquad,&nquad,NULL);CHKERRQ(ierr);
   ierr        = PetscOptionsEnd();
   user.delxyz  = 1.0/(user.n);
+  ierr = PetscMalloc(sizeof(PetscReal)*nquad, &user.quadx);CHKERRQ(ierr);
+  ierr = PetscMalloc(sizeof(PetscReal)*nquad, &user.quadw);CHKERRQ(ierr);
 
   /* create the DM object for dimension specified */
-  ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, NULL, user.n, 1, &dm);CHKERRQ(ierr);
+  const double bounds[4] = {0.0,2.5,0.0,2.5};
+  ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, bounds, user.n, 1, &dm);CHKERRQ(ierr);
   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
   ierr = DMMoabSetBlockSize(dm, 2);CHKERRQ(ierr);
   ierr = DMMoabSetFields(dm, 2, fields);CHKERRQ(ierr);
   /* SetUp the data structures for DMMOAB */
   ierr = DMSetUp(dm);CHKERRQ(ierr);
 
+  /* Compute the Gauss quadrature points and corresponding weights */
+  ierr = PetscDTGaussQuadrature(nquad, -1.0, 1.0, user.quadx, user.quadw);CHKERRQ(ierr);
+
   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      Extract global vectors from DMDA; then duplicate for remaining
      vectors that are the same types
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   ierr = DMCreateGlobalVector(dm,&x);CHKERRQ(ierr);
-//  ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
   ierr = DMCreateMatrix(dm,MATBAIJ,&A);CHKERRQ(ierr);
 
   mon.comm    = PETSC_COMM_WORLD;
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
   ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
-  ierr = TSSetType(ts,TSGL);CHKERRQ(ierr); /* Rosenbrock-W */
+  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr); /* Rosenbrock-W */
   ierr = TSSetIFunction(ts,NULL,NeutHC_IFunction,&user);CHKERRQ(ierr);
   ierr = TSSetIJacobian(ts,A,A,NeutHC_IJacobian,&user);CHKERRQ(ierr);
-
   
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      Create matrix data structure; set Jacobian evaluation routine
      -snes_mf_operator : form preconditioning matrix as set by the user,
                          but use matrix-free approx for Jacobian-vector
                          products within Newton-Krylov method
-
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
   ierr = TSSetDuration(ts,maxsteps,endtime);CHKERRQ(ierr);
      Set initial conditions
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   ierr = NeutHC_InitialSolution(dm,0.0,x,&user);CHKERRQ(ierr);
-  ierr = TSSetInitialTimeStep(ts,0.0,.0001);CHKERRQ(ierr);
+  ierr = TSSetInitialTimeStep(ts,0.0,1.0e-8);CHKERRQ(ierr);
   ierr = TSSetSolution(ts,x);CHKERRQ(ierr);
 
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   PetscScalar       gradU,evolume;
   PetscInt          i,j,q,num_conn,offset;
   const moab::EntityHandle *connect;
-  moab::Range       elocal;
+  const moab::Range *elocal;
   moab::Interface*  mbImpl;
   Field             locala[VPERE*VPERE],localr[VPERE];
   PetscBool         elem_on_boundary;
   ierr = DMMoabGetLocalToGlobalTag(dm, &idtag);CHKERRQ(ierr); 
 
   /* loop over local elements */
-  for(moab::Range::iterator iter = elocal.begin(); iter != elocal.end(); iter++) {
+  for(moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
     const moab::EntityHandle ehandle = *iter;
 
     ierr = PetscMemzero(localr,sizeof(Field)*VPERE);CHKERRQ(ierr);
 
     merr = mbImpl->tag_get_data(idtag,&ehandle,1,&egid);MBERRNM(merr);
 
-//    PetscPrintf(PETSC_COMM_WORLD, "ELEMENT %D \t TOTALXS = %G \t AVGTEMP = %G\n",egid,totalxs,avgtemp);
+//    PetscPrintf(PETSC_COMM_WORLD, "ELEMENT %D \t DOFS = [%D, %D, %D, %D]\n",egid,dof_indices[0], dof_indices[1], dof_indices[2], dof_indices[3]);
 
     /* Compute function over the locally owned part of the grid */
     for (q=0; q<NQPTS; ++q) {
       gradU=0.0;
       xqpts.U=xqpts.V=0.0;
       xdotqpts.U=xdotqpts.V=0.0;
+      offset=q*VPERE;
       for (i=0; i < VPERE; ++i) {
-        xqpts.U+=phi[q*VPERE+i]*x[dof_indices[i]].U;
-        xqpts.V+=phi[q*VPERE+i]*x[dof_indices[i]].V;
-        xdotqpts.U+= phi[q*VPERE+i]*xdot[dof_indices[i]].U;
-        xdotqpts.V+= phi[q*VPERE+i]*xdot[dof_indices[i]].V;
+        xqpts.U+=phi[offset+i]*x[dof_indices[i]].U;
+        xqpts.V+=phi[offset+i]*x[dof_indices[i]].V;
+        xdotqpts.U+= phi[offset+i]*xdot[dof_indices[i]].U;
+        xdotqpts.V+= phi[offset+i]*xdot[dof_indices[i]].V;
       }
+//    PetscPrintf(PETSC_COMM_WORLD, "\t Quadrature %D \t [U, V] = [%G, %G, %G, %G]\n",egid,xqpts.U, xqpts.V, xdotqpts.U, xdotqpts.V);
 
       for (i=0; i < VPERE; ++i) {
+        offset=i*VPERE;
         for (j=0; j < VPERE; ++j) {
           gradU += ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) ;
 
-          locala[i*VPERE+j].U += jxw[q] * ( user->diffusion.U * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) +
-                                              (user->eta + xqpts.V * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
+          locala[offset+j].U += jxw[q] * ( user->diffusion.U * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) +
+                                              (xdotqpts.U + user->eta + xqpts.V * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
 
-          locala[i*VPERE+j].V += jxw[q] * ( user->diffusion.V * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) + 
-                                              (user->eta + user->psi - xqpts.U * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
+          locala[offset+j].V += jxw[q] * ( user->diffusion.V * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) + 
+                                              (xdotqpts.V + user->eta + user->psi - xqpts.U * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
 
         }
-//        localf[i].U += jxw[q] * phi[q*VPERE+i] * (-xdot[dof_indices[i]].U);
+
 #if (SUBPROBLEM==1)
-        localr[i].U += jxw[q] * phi[q*VPERE+i] * (xdotqpts.U - user->eta);
+        localr[i].U += jxw[q] * phi[q*VPERE+i] * ( - user->eta);
 #else
         localr[i].U += jxw[q] * phi[q*VPERE+i] * (xdotqpts.U - user->eta + user->cpld_diffusion * gradU * x[dof_indices[i]].U);
 #endif
-        localr[i].V += jxw[q] * phi[q*VPERE+i] * xdotqpts.V ;
+//        localr[i].V += jxw[q] * phi[q*VPERE+i] * (xdotqpts.V);
       }
     }
 
       array[dof_indices[i]].U += localr[i].U;
       array[dof_indices[i]].V += localr[i].V;
     }
+//  ierr = DMMoabVecRestoreArray(dm, F, &array);CHKERRQ(ierr);
+//  VecView(F,0);
+//  ierr = DMMoabVecGetArray(dm, F, &array);CHKERRQ(ierr);    
 #else
     ierr = VecSetValuesBlocked(F, VPERE, dof_indices, (PetscScalar*)localr, ADD_VALUES);CHKERRQ(ierr);
 #endif
   PetscScalar       evolume;
   PetscInt          i,j,q,num_conn;
   const moab::EntityHandle *connect;
-  moab::Range       elocal;
+  const moab::Range *elocal;
   moab::Interface*  mbImpl;
-  Field             locala[VPERE*VPERE];
+  PetscScalar       locala[VPERE][2][VPERE][2];
   PetscErrorCode    ierr;
   PetscScalar       phi[VPERE*NQPTS],dphidx[VPERE*NQPTS],dphidy[VPERE*NQPTS],dphidz[VPERE*NQPTS];
   PetscScalar       vpos[VPERE*3],quadrature[NQPTS*3],jxw[NQPTS];
   ierr = DMMoabGetLocalToGlobalTag(dm, &idtag);CHKERRQ(ierr);
 
   /* loop over local elements */
-  for(moab::Range::iterator iter = elocal.begin(); iter != elocal.end(); iter++) {
+  for(moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
     const moab::EntityHandle ehandle = *iter;
 
-    ierr = PetscMemzero(locala,sizeof(Field)*VPERE*VPERE);CHKERRQ(ierr);
+    ierr = PetscMemzero(locala,sizeof(locala));CHKERRQ(ierr);
 
     // Get connectivity information:
     ierr = DMMoabGetElementConnectivity(dm, ehandle, &num_conn, &connect);CHKERRQ(ierr);
     ierr = DMMoabGetVertexCoordinates(dm,num_conn,connect,vpos);CHKERRQ(ierr);
 
     /* get the global DOF number to appropriately set the element contribution in the RHS vector */
-    ierr = DMMoabGetDofsBlocked(dm, num_conn, connect, dof_indices);CHKERRQ(ierr);
+    ierr = DMMoabGetDofsBlockedLocal(dm, num_conn, connect, dof_indices);CHKERRQ(ierr);
 
     /* compute the quadrature points transformed to the physical space */
     ierr = ComputeQuadraturePointsPhysical(user->dim, vpos, quadrature, jxw);CHKERRQ(ierr);
 
       for (i=0; i < VPERE; ++i) {
         for (j=0; j < VPERE; ++j) {
-          locala[i*VPERE+j].U += jxw[q] * ( user->diffusion.U * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) +
-                                              (afactor + user->eta + xqpts.V * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
+          locala[i][0][j][0] += jxw[q] * ( user->diffusion.U * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) +
+                                              (afactor + user->eta) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
+
+          locala[i][0][j][1] += jxw[q] * ( (xqpts.V * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
 
-          locala[i*VPERE+j].V += jxw[q] * ( user->diffusion.V * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) + 
-                                              (afactor + user->eta + user->psi - xqpts.U * xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
+          locala[i][1][j][1] += jxw[q] * ( user->diffusion.V * ( dphidx[q*VPERE+i]*dphidx[q*VPERE+j] + dphidy[q*VPERE+i]*dphidy[q*VPERE+j] ) + 
+                                              (afactor + user->eta + user->psi - xqpts.V) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
+
+          locala[i][1][j][0] += jxw[q] * ( (- xqpts.U) * phi[q*VPERE+i] * phi[q*VPERE+j] ) ;
         }
       }
     }
 
     /* set the values directly into appropriate locations. Can alternately use VecSetValues */
-    ierr = MatSetValuesBlocked(*B, VPERE, dof_indices, VPERE, dof_indices, (PetscScalar*)locala, ADD_VALUES);CHKERRQ(ierr);
+    ierr = MatSetValuesBlockedLocal(*B, VPERE, dof_indices, VPERE, dof_indices, &locala[0][0][0][0], ADD_VALUES);CHKERRQ(ierr);
   }
 
   ierr = DMMoabVecRestoreArrayRead(dm, X, &x);CHKERRQ(ierr);
   Field             *x;
   PetscErrorCode    ierr;
   moab::Interface*  mbImpl;
-  moab::Range       vowned;
+  const moab::Range *vowned;
   PetscInt          dof_index;
   moab::Range::iterator iter;
   PetscReal         vpos[3];
   PetscRandom       rnd;
+  const PetscReal   h=1.25,k=1.25;
   PetscScalar       value;
 #if (SUBPROBLEM==1)
   const PetscScalar u0 = 1.0;
 #elif defined(PETSC_HAVE_RAND)
   ierr = PetscRandomSetType(rnd,PETSCRAND);CHKERRQ(ierr);
 #endif
-  ierr = PetscRandomSetInterval(rnd, -0.01, 0.01);CHKERRQ(ierr);
+  ierr = PetscRandomSetInterval(rnd, -0.05, 0.05);CHKERRQ(ierr);
   ierr = PetscRandomSetFromOptions(rnd);CHKERRQ(ierr);
 
   /* Compute function over the locally owned part of the grid */
-  for(moab::Range::iterator iter = vowned.begin(); iter != vowned.end(); iter++) {
+  for(moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) {
     const moab::EntityHandle vhandle = *iter;
     ierr = DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index);CHKERRQ(ierr);
 
     /* compute the mid-point of the element and use a 1-point lumped quadrature */
     ierr = DMMoabGetVertexCoordinates(dm,1,&vhandle,vpos);CHKERRQ(ierr);
 
-    if ((vpos[0] >= 0.45 && vpos[0] <= 0.55) && (vpos[1] >= 0.45 && vpos[1] <= 0.55)) {
+    ierr = PetscRandomGetValue(rnd, &value);CHKERRQ(ierr);
+#if (SUBPROBLEM!=1)
+    x[dof_index].U = u0 - 1.0/10.0 * PetscExpScalar(-50.0*((vpos[0]-h)*(vpos[0]-h)+(vpos[1]-k)*(vpos[1]-k)));
+    x[dof_index].V = v0 + 1.0/10.0 * PetscExpScalar(-50.0*((vpos[0]-h)*(vpos[0]-h)+(vpos[1]-k)*(vpos[1]-k)));
+#else
       /* Reference: Pearson, "Complex Patterns in a Simple System", http://arxiv.org/pdf/patt-sol/9304003.pdf */
-      ierr = PetscRandomGetValue(rnd, &value);CHKERRQ(ierr);
+    if ((vpos[0] >= 1.0 && vpos[0] <= 1.5) && (vpos[1] >= 1.0 && vpos[1] <= 1.5)) {
 //      x[dof_index].U = u0 * (1.0 + (value/100000.0));
-      x[dof_index].U = 0.5 * (1.0 + value);;
-      x[dof_index].V = 0.25 * (1.0 + value);;
+      x[dof_index].U = 0.5 * (1.0 + value) * PetscExpScalar(-20.0*((vpos[0]-h)*(vpos[0]-h)+(vpos[1]-k)*(vpos[1]-k)));
+      x[dof_index].V = 0.25 * (1.0 + value) * PetscExpScalar(-20.0*((vpos[0]-h)*(vpos[0]-h)+(vpos[1]-k)*(vpos[1]-k)));
     }
     else {
       x[dof_index].U = u0;
       x[dof_index].V = v0;
     }
+#endif
   }
 
   ierr = PetscRandomDestroy(&rnd);CHKERRQ(ierr);
         dphidx[0+j*2] = -1.0;
         dphidx[1+j*2] = 1.0;
       }
-
-      if (dphidy) {
-        dphidy[0+j*2] = dphidy[1+j*2] = 0.0;
-      }
     }
   }
   else if (dim==2) {
       phi[i+j*VPERE]    *= ejac;
       if (dphidx) dphidx[i+j*VPERE] *= ejac;
       if (dphidy) dphidy[i+j*VPERE] *= ejac;
+      if (dphidz) dphidz[i+j*VPERE] *= ejac;
     }
   }
 
 #if 0
   /* verify if the computed basis functions are consistent */
   for ( j = 0; j < n; j++ ) {
-    PetscScalar phisum=0,dphixsum=0,dphiysum=0;
+    PetscScalar phisum=0,dphixsum=0,dphiysum=0,dphizsum=0;
     for ( i = 0; i < VPERE; i++ ) {
       phisum += phi[i+j*4];
       if (dphidx) dphixsum += dphidx[i+j*4];
       if (dphidy) dphiysum += dphidy[i+j*4];
+      if (dphidz) dphizsum += dphidz[i+j*4];
     }
-    PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %G, %G, %G\n", j, phisum, dphixsum, dphiysum);
+    PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %G, %G, %G\n", j, phisum, dphixsum, dphiysum, dphizsum);
   }
 #endif
   PetscFunctionReturn(0);