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>
5
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;
17
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;
23
24 /*
25 FormPermeability - Forms permeability field.
26
27 Input Parameters:
28 user - user-defined application context
29
30 Output Parameter:
31 Kappa - vector
32 */
FormPermeability(DM da,Vec Kappa,AppCtx * user)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;
40
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 }
60
61 /*
62 FormFunctionLocal - Evaluates nonlinear residual, F(x) on local process patch
63 */
FormFunctionLocal(DMDALocalInfo * info,Field * u,Field * f,AppCtx * user)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;
77
78 PetscFunctionBeginUser;
79 PetscCall(DMGetGlobalVector(user->cda, &L));
80
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));
95
96 f[i].s = phi * (u[i].s - uold[i].s) + (dt / dx) * ((lambdaWet / lambda) * u[i].v - (lambdaWetL / lambdaL) * u[i - 1].v);
97
98 f[i].v = u[i].v + K * lambda * (u[i].p - u[i - 1].p) / dx;
99
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)); */
107
108 PetscCall(DMRestoreGlobalVector(user->cda, &L));
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
main(int argc,char ** argv)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;
120
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));
142
143 user.sl = 1.0;
144 user.vl = 0.1;
145 user.pl = 1.0;
146 user.phi = 1.0;
147
148 user.kappaWet = 1.0;
149 user.kappaNoWet = 0.3;
150
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));
161
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 }
174
175 /*TEST
176
177 test:
178 suffix: 0
179 requires: !single
180 args: -snes_converged_reason -snes_monitor_short
181
182 TEST*/
183