Actual source code: ex33.c
petsc-3.9.4 2018-09-11
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;
43: DMGetCoordinateDM(da, &cda);
44: DMGetCoordinates(da, &c);
45: DMDAGetCorners(da, &xs,NULL,NULL, &xm,NULL,NULL);
46: DMDAVecGetArray(da, Kappa, &K);
47: DMDAVecGetArray(cda, c, &coords);
48: for (i = xs; i < xs+xm; ++i) {
49: #if 1
50: K[i] = 1.0;
51: #else
52: /* Notch */
53: if (i == (xs+xm)/2) K[i] = 0.00000001;
54: else K[i] = 1.0;
55: #endif
56: }
57: DMDAVecRestoreArray(da, Kappa, &K);
58: DMDAVecRestoreArray(cda, c, &coords);
59: return(0);
60: }
62: /*
63: FormFunctionLocal - Evaluates nonlinear residual, F(x) on local process patch
64: */
65: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field *u, Field *f, AppCtx *user)
66: {
67: Vec L;
68: PetscReal phi = user->phi;
69: PetscReal dt = user->dt;
70: PetscReal dx = 1.0/(PetscReal)(info->mx-1);
71: PetscReal alpha = 2.0;
72: PetscReal beta = 2.0;
73: PetscReal kappaWet = user->kappaWet;
74: PetscReal kappaNoWet = user->kappaNoWet;
75: Field *uold;
76: PetscScalar *Kappa;
77: PetscInt i;
81: DMGetGlobalVector(user->cda, &L);
83: DMDAVecGetArray(info->da, user->uold, &uold);
84: DMDAVecGetArray(user->cda, user->Kappa, &Kappa);
85: /* Compute residual over the locally owned part of the grid */
86: for (i = info->xs; i < info->xs+info->xm; ++i) {
87: if (i == 0) {
88: f[i].s = u[i].s - user->sl;
89: f[i].v = u[i].v - user->vl;
90: f[i].p = u[i].p - user->pl;
91: } else {
92: PetscScalar K = 2*dx/(dx/Kappa[i] + dx/Kappa[i-1]);
93: PetscReal lambdaWet = kappaWet*PetscRealPart(PetscPowScalar(u[i].s, alpha));
94: PetscReal lambda = lambdaWet + kappaNoWet*PetscRealPart(PetscPowScalar(1.-u[i].s, beta));
95: PetscReal lambdaWetL = kappaWet*PetscRealPart(PetscPowScalar(u[i-1].s, alpha));
96: PetscReal lambdaL = lambdaWetL + kappaNoWet*PetscRealPart(PetscPowScalar(1.-u[i-1].s, beta));
98: f[i].s = phi*(u[i].s - uold[i].s) + (dt/dx)*((lambdaWet/lambda)*u[i].v - (lambdaWetL/lambdaL)*u[i-1].v);
100: f[i].v = u[i].v + K*lambda*(u[i].p - u[i-1].p)/dx;
102: /*pxx = (2.0*u[i].p - u[i-1].p - u[i+1].p)/dx;*/
103: f[i].p = u[i].v - u[i-1].v;
104: }
105: }
106: DMDAVecRestoreArray(info->da, user->uold, &uold);
107: DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa);
108: /* PetscLogFlops(11.0*info->ym*info->xm); */
110: DMRestoreGlobalVector(user->cda, &L);
111: return(0);
112: }
114: int main(int argc, char **argv)
115: {
116: SNES snes; /* nonlinear solver */
117: DM da; /* grid */
118: Vec u; /* solution vector */
119: AppCtx user; /* user-defined work context */
120: PetscReal t = 0.0;/* time */
122: PetscInt n;
124: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
125: /* Create solver */
126: SNESCreate(PETSC_COMM_WORLD, &snes);
127: /* Create mesh */
128: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,3,1,NULL,&da);
129: DMSetFromOptions(da);
130: DMSetUp(da);
131: DMSetApplicationContext(da, &user);
132: SNESSetDM(snes, da);
133: /* Create coefficient */
134: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,1,1,NULL,&user.cda);
135: DMSetFromOptions(user.cda);
136: DMSetUp(user.cda);
137: DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
138: DMGetGlobalVector(user.cda, &user.Kappa);
139: FormPermeability(user.cda, user.Kappa, &user);
140: /* Setup Problem */
141: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
142: DMGetGlobalVector(da, &u);
143: DMGetGlobalVector(da, &user.uold);
145: user.sl = 1.0;
146: user.vl = 0.1;
147: user.pl = 1.0;
148: user.phi = 1.0;
150: user.kappaWet = 1.0;
151: user.kappaNoWet = 0.3;
153: /* Time Loop */
154: user.dt = 0.1;
155: for (n = 0; n < 100; ++n, t += user.dt) {
156: PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", (double)t);
157: VecView(u, PETSC_VIEWER_DRAW_WORLD);
158: /* Solve */
159: SNESSetFromOptions(snes);
160: SNESSolve(snes, NULL, u);
161: /* Update */
162: VecCopy(u, user.uold);
164: VecView(u, PETSC_VIEWER_DRAW_WORLD);
165: }
166: /* Cleanup */
167: DMRestoreGlobalVector(da, &u);
168: DMRestoreGlobalVector(da, &user.uold);
169: DMRestoreGlobalVector(user.cda, &user.Kappa);
170: DMDestroy(&user.cda);
171: DMDestroy(&da);
172: SNESDestroy(&snes);
173: PetscFinalize();
174: return ierr;
175: }
177: /*TEST
179: test:
180: suffix: 0
181: requires: !single
182: args: -snes_converged_reason -snes_monitor_short
184: TEST*/