1: static char help[] = "Multiphase flow in a porous medium in 1d.\n\n";
  2: #include <petscdm.h>
  3: #include <petscdmda.h>
  4: #include <petscsnes.h>

  6: typedef struct {
  7:   DM        cda;
  8:   Vec       uold;
  9:   Vec       Kappa;
 10:   PetscReal phi;
 11:   PetscReal kappaWet;
 12:   PetscReal kappaNoWet;
 13:   PetscReal dt;
 14:   /* Boundary conditions */
 15:   PetscReal sl, vl, pl;
 16: } AppCtx;

 18: typedef struct {
 19:   PetscScalar s; /* The saturation on each cell */
 20:   PetscScalar v; /* The velocity on each face */
 21:   PetscScalar p; /* The pressure on each cell */
 22: } Field;

 26: /*
 27:    FormPermeability - Forms permeability field.

 29:    Input Parameters:
 30:    user - user-defined application context

 32:    Output Parameter:
 33:    Kappa - vector
 34:  */
 35: PetscErrorCode FormPermeability(DM da, Vec Kappa, AppCtx *user)
 36: {
 37:   DM             cda;
 38:   Vec            c;
 39:   PetscScalar    *K;
 40:   PetscScalar    *coords;
 41:   PetscInt       xs, xm, i;

 45:   DMGetCoordinateDM(da, &cda);
 46:   DMGetCoordinates(da, &c);
 47:   DMDAGetCorners(da, &xs,NULL,NULL, &xm,NULL,NULL);
 48:   DMDAVecGetArray(da, Kappa, &K);
 49:   DMDAVecGetArray(cda, c, &coords);
 50:   for (i = xs; i < xs+xm; ++i) {
 51: #if 1
 52:     K[i] = 1.0;
 53: #else
 54:     /* Notch */
 55:     if (i == (xs+xm)/2) K[i] = 0.00000001;
 56:     else K[i] = 1.0;
 57: #endif
 58:   }
 59:   DMDAVecRestoreArray(da, Kappa, &K);
 60:   DMDAVecRestoreArray(cda, c, &coords);
 61:   return(0);
 62: }

 66: /*
 67:    FormFunctionLocal - Evaluates nonlinear residual, F(x) on local process patch
 68: */
 69: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field *u, Field *f, AppCtx *user)
 70: {
 71:   Vec            L;
 72:   PetscReal      phi        = user->phi;
 73:   PetscReal      dt         = user->dt;
 74:   PetscReal      dx         = 1.0/(PetscReal)(info->mx-1);
 75:   PetscReal      alpha      = 2.0;
 76:   PetscReal      beta       = 2.0;
 77:   PetscReal      kappaWet   = user->kappaWet;
 78:   PetscReal      kappaNoWet = user->kappaNoWet;
 79:   Field          *uold;
 80:   PetscScalar    *Kappa;
 81:   PetscInt       i;

 85:   DMGetGlobalVector(user->cda, &L);

 87:   DMDAVecGetArray(info->da, user->uold,  &uold);
 88:   DMDAVecGetArray(user->cda, user->Kappa, &Kappa);
 89:   /* Compute residual over the locally owned part of the grid */
 90:   for (i = info->xs; i < info->xs+info->xm; ++i) {
 91:     if (i == 0) {
 92:       f[i].s = u[i].s - user->sl;
 93:       f[i].v = u[i].v - user->vl;
 94:       f[i].p = u[i].p - user->pl;
 95:     } else {
 96:       PetscScalar K          = 2*dx/(dx/Kappa[i] + dx/Kappa[i-1]);
 97:       PetscReal   lambdaWet  = kappaWet*PetscPowScalar(u[i].s, alpha);
 98:       PetscReal   lambda     = lambdaWet + kappaNoWet*PetscPowScalar(1-u[i].s, beta);
 99:       PetscReal   lambdaWetL = kappaWet*PetscPowScalar(u[i-1].s, alpha);
100:       PetscReal   lambdaL    = lambdaWetL + kappaNoWet*PetscPowScalar(1-u[i-1].s, beta);

102:       f[i].s = phi*(u[i].s - uold[i].s) + (dt/dx)*((lambdaWet/lambda)*u[i].v - (lambdaWetL/lambdaL)*u[i-1].v);

104:       f[i].v = u[i].v + K*lambda*(u[i].p - u[i-1].p)/dx;

106:       /*pxx     = (2.0*u[i].p - u[i-1].p - u[i+1].p)/dx;*/
107:       f[i].p = u[i].v - u[i-1].v;
108:     }
109:   }
110:   DMDAVecRestoreArray(info->da, user->uold, &uold);
111:   DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa);
112:   /* PetscLogFlops(11.0*info->ym*info->xm); */

114:   DMRestoreGlobalVector(user->cda, &L);
115:   return(0);
116: }

120: int main(int argc, char **argv)
121: {
122:   SNES           snes;   /* nonlinear solver */
123:   DM             da;     /* grid */
124:   Vec            u;      /* solution vector */
125:   AppCtx         user;   /* user-defined work context */
126:   PetscReal      t;      /* time */
128:   PetscInt       n;

130:   PetscInitialize(&argc, &argv, NULL, help);
131:   /* Create solver */
132:   SNESCreate(PETSC_COMM_WORLD, &snes);
133:   /* Create mesh */
134:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-4,3,1,NULL,&da);
135:   DMSetApplicationContext(da, &user);
136:   SNESSetDM(snes, da);
137:   /* Create coefficient */
138:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-4,1,1,NULL,&user.cda);
139:   DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
140:   DMGetGlobalVector(user.cda, &user.Kappa);
141:   FormPermeability(user.cda, user.Kappa, &user);
142:   /* Setup Problem */
143:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
144:   DMGetGlobalVector(da, &u);
145:   DMGetGlobalVector(da, &user.uold);

147:   user.sl  = 1.0;
148:   user.vl  = 0.1;
149:   user.pl  = 1.0;
150:   user.phi = 1.0;

152:   user.kappaWet   = 1.0;
153:   user.kappaNoWet = 0.3;

155:   /* Time Loop */
156:   user.dt = 0.1;
157:   for (n = 0; n < 100; ++n, t += user.dt) {
158:     PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", t);
159:     VecView(u, PETSC_VIEWER_DRAW_WORLD);
160:     /* Solve */
161:     SNESSetFromOptions(snes);
162:     SNESSolve(snes, NULL, u);
163:     /* Update */
164:     VecCopy(u, user.uold);

166:     VecView(u, PETSC_VIEWER_DRAW_WORLD);
167:   }
168:   /* Cleanup */
169:   DMRestoreGlobalVector(da, &u);
170:   DMRestoreGlobalVector(da, &user.uold);
171:   DMRestoreGlobalVector(user.cda, &user.Kappa);
172:   DMDestroy(&user.cda);
173:   DMDestroy(&da);
174:   SNESDestroy(&snes);
175:   PetscFinalize();
176:   return 0;
177: }