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;
77: PetscFunctionBeginUser;
78: PetscCall(DMGetGlobalVector(user->cda, &L));
80: PetscCall(DMDAVecGetArray(info->da, user->uold, &uold));
81: PetscCall(DMDAVecGetArray(user->cda, user->Kappa, &Kappa));
82: /* Compute residual over the locally owned part of the grid */
83: for (PetscInt i = info->xs; i < info->xs + info->xm; ++i) {
84: if (i == 0) {
85: f[i].s = u[i].s - user->sl;
86: f[i].v = u[i].v - user->vl;
87: f[i].p = u[i].p - user->pl;
88: } else {
89: PetscScalar K = 2 * dx / (dx / Kappa[i] + dx / Kappa[i - 1]);
90: PetscReal lambdaWet = kappaWet * PetscRealPart(PetscPowScalar(u[i].s, alpha));
91: PetscReal lambda = lambdaWet + kappaNoWet * PetscRealPart(PetscPowScalar(1. - u[i].s, beta));
92: PetscReal lambdaWetL = kappaWet * PetscRealPart(PetscPowScalar(u[i - 1].s, alpha));
93: PetscReal lambdaL = lambdaWetL + kappaNoWet * PetscRealPart(PetscPowScalar(1. - u[i - 1].s, beta));
95: f[i].s = phi * (u[i].s - uold[i].s) + (dt / dx) * ((lambdaWet / lambda) * u[i].v - (lambdaWetL / lambdaL) * u[i - 1].v);
97: f[i].v = u[i].v + K * lambda * (u[i].p - u[i - 1].p) / dx;
99: /*pxx = (2.0*u[i].p - u[i-1].p - u[i+1].p)/dx;*/
100: f[i].p = u[i].v - u[i - 1].v;
101: }
102: }
103: PetscCall(DMDAVecRestoreArray(info->da, user->uold, &uold));
104: PetscCall(DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa));
105: /* PetscCall(PetscLogFlops(11.0*info->ym*info->xm)); */
107: PetscCall(DMRestoreGlobalVector(user->cda, &L));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: int main(int argc, char **argv)
112: {
113: SNES snes; /* nonlinear solver */
114: DM da; /* grid */
115: Vec u; /* solution vector */
116: AppCtx user; /* user-defined work context */
117: PetscReal t = 0.0; /* time */
119: PetscFunctionBeginUser;
120: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
121: /* Create solver */
122: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
123: /* Create mesh */
124: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 3, 1, NULL, &da));
125: PetscCall(DMSetFromOptions(da));
126: PetscCall(DMSetUp(da));
127: PetscCall(DMSetApplicationContext(da, &user));
128: PetscCall(SNESSetDM(snes, da));
129: /* Create coefficient */
130: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 1, 1, NULL, &user.cda));
131: PetscCall(DMSetFromOptions(user.cda));
132: PetscCall(DMSetUp(user.cda));
133: PetscCall(DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
134: PetscCall(DMGetGlobalVector(user.cda, &user.Kappa));
135: PetscCall(FormPermeability(user.cda, user.Kappa, &user));
136: /* Setup Problem */
137: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, &user));
138: PetscCall(DMGetGlobalVector(da, &u));
139: PetscCall(DMGetGlobalVector(da, &user.uold));
141: user.sl = 1.0;
142: user.vl = 0.1;
143: user.pl = 1.0;
144: user.phi = 1.0;
146: user.kappaWet = 1.0;
147: user.kappaNoWet = 0.3;
149: /* Time Loop */
150: user.dt = 0.1;
151: for (PetscInt n = 0; n < 100; ++n, t += user.dt) {
152: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", (double)t));
153: PetscCall(VecView(u, PETSC_VIEWER_DRAW_WORLD));
154: /* Solve */
155: PetscCall(SNESSetFromOptions(snes));
156: PetscCall(SNESSolve(snes, NULL, u));
157: /* Update */
158: PetscCall(VecCopy(u, user.uold));
160: PetscCall(VecView(u, PETSC_VIEWER_DRAW_WORLD));
161: }
162: /* Cleanup */
163: PetscCall(DMRestoreGlobalVector(da, &u));
164: PetscCall(DMRestoreGlobalVector(da, &user.uold));
165: PetscCall(DMRestoreGlobalVector(user.cda, &user.Kappa));
166: PetscCall(DMDestroy(&user.cda));
167: PetscCall(DMDestroy(&da));
168: PetscCall(SNESDestroy(&snes));
169: PetscCall(PetscFinalize());
170: return 0;
171: }
173: /*TEST
175: test:
176: suffix: 0
177: requires: !single
178: args: -snes_converged_reason -snes_monitor_short
180: TEST*/