Actual source code: ex33.c

      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;

     24: /*
     25:    FormPermeability - Forms permeability field.

     27:    Input Parameters:
     28:    user - user-defined application context

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

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

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

     78:   PetscFunctionBeginUser;
     79:   PetscCall(DMGetGlobalVector(user->cda, &L));

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

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

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

    100:       /*pxx     = (2.0*u[i].p - u[i-1].p - u[i+1].p)/dx;*/
    101:       f[i].p = u[i].v - u[i - 1].v;
    102:     }
    103:   }
    104:   PetscCall(DMDAVecRestoreArray(info->da, user->uold, &uold));
    105:   PetscCall(DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa));
    106:   /* PetscCall(PetscLogFlops(11.0*info->ym*info->xm)); */

    108:   PetscCall(DMRestoreGlobalVector(user->cda, &L));
    109:   PetscFunctionReturn(PETSC_SUCCESS);
    110: }

    112: int main(int argc, char **argv)
    113: {
    114:   SNES      snes;    /* nonlinear solver */
    115:   DM        da;      /* grid */
    116:   Vec       u;       /* solution vector */
    117:   AppCtx    user;    /* user-defined work context */
    118:   PetscReal t = 0.0; /* time */
    119:   PetscInt  n;

    121:   PetscFunctionBeginUser;
    122:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    123:   /* Create solver */
    124:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    125:   /* Create mesh */
    126:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 3, 1, NULL, &da));
    127:   PetscCall(DMSetFromOptions(da));
    128:   PetscCall(DMSetUp(da));
    129:   PetscCall(DMSetApplicationContext(da, &user));
    130:   PetscCall(SNESSetDM(snes, da));
    131:   /* Create coefficient */
    132:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 1, 1, NULL, &user.cda));
    133:   PetscCall(DMSetFromOptions(user.cda));
    134:   PetscCall(DMSetUp(user.cda));
    135:   PetscCall(DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
    136:   PetscCall(DMGetGlobalVector(user.cda, &user.Kappa));
    137:   PetscCall(FormPermeability(user.cda, user.Kappa, &user));
    138:   /* Setup Problem */
    139:   PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
    140:   PetscCall(DMGetGlobalVector(da, &u));
    141:   PetscCall(DMGetGlobalVector(da, &user.uold));

    143:   user.sl  = 1.0;
    144:   user.vl  = 0.1;
    145:   user.pl  = 1.0;
    146:   user.phi = 1.0;

    148:   user.kappaWet   = 1.0;
    149:   user.kappaNoWet = 0.3;

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

    162:     PetscCall(VecView(u, PETSC_VIEWER_DRAW_WORLD));
    163:   }
    164:   /* Cleanup */
    165:   PetscCall(DMRestoreGlobalVector(da, &u));
    166:   PetscCall(DMRestoreGlobalVector(da, &user.uold));
    167:   PetscCall(DMRestoreGlobalVector(user.cda, &user.Kappa));
    168:   PetscCall(DMDestroy(&user.cda));
    169:   PetscCall(DMDestroy(&da));
    170:   PetscCall(SNESDestroy(&snes));
    171:   PetscCall(PetscFinalize());
    172:   return 0;
    173: }

    175: /*TEST

    177:   test:
    178:     suffix: 0
    179:     requires: !single
    180:     args: -snes_converged_reason -snes_monitor_short

    182: TEST*/