Actual source code: ex33.c
petsc-3.4.5 2014-06-29
1: static char help[] = "Multiphase flow in a porous medium in 1d.\n\n";
2: #include <petscdmda.h>
3: #include <petscsnes.h>
5: typedef struct {
6: DM cda;
7: Vec uold;
8: Vec Kappa;
9: PetscReal phi;
10: PetscReal kappaWet;
11: PetscReal kappaNoWet;
12: PetscReal dt;
13: /* Boundary conditions */
14: PetscReal sl, vl, pl;
15: } AppCtx;
17: typedef struct {
18: PetscScalar s; /* The saturation on each cell */
19: PetscScalar v; /* The velocity on each face */
20: PetscScalar p; /* The pressure on each cell */
21: } Field;
25: /*
26: FormPermeability - Forms permeability field.
28: Input Parameters:
29: user - user-defined application context
31: Output Parameter:
32: Kappa - vector
33: */
34: PetscErrorCode FormPermeability(DM da, Vec Kappa, AppCtx *user)
35: {
36: DM cda;
37: Vec c;
38: PetscScalar *K;
39: PetscScalar *coords;
40: PetscInt xs, xm, i;
44: DMGetCoordinateDM(da, &cda);
45: DMGetCoordinates(da, &c);
46: DMDAGetCorners(da, &xs,NULL,NULL, &xm,NULL,NULL);
47: DMDAVecGetArray(da, Kappa, &K);
48: DMDAVecGetArray(cda, c, &coords);
49: for (i = xs; i < xs+xm; ++i) {
50: #if 1
51: K[i] = 1.0;
52: #else
53: /* Notch */
54: if (i == (xs+xm)/2) K[i] = 0.00000001;
55: else K[i] = 1.0;
56: #endif
57: }
58: DMDAVecRestoreArray(da, Kappa, &K);
59: DMDAVecRestoreArray(cda, c, &coords);
60: return(0);
61: }
65: /*
66: FormFunctionLocal - Evaluates nonlinear residual, F(x) on local process patch
67: */
68: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field *u, Field *f, AppCtx *user)
69: {
70: Vec L;
71: PetscReal phi = user->phi;
72: PetscReal dt = user->dt;
73: PetscReal dx = 1.0/(PetscReal)(info->mx-1);
74: PetscReal alpha = 2.0;
75: PetscReal beta = 2.0;
76: PetscReal kappaWet = user->kappaWet;
77: PetscReal kappaNoWet = user->kappaNoWet;
78: Field *uold;
79: PetscScalar *Kappa;
80: PetscInt i;
84: DMGetGlobalVector(user->cda, &L);
86: DMDAVecGetArray(info->da, user->uold, &uold);
87: DMDAVecGetArray(user->cda, user->Kappa, &Kappa);
88: /* Compute residual over the locally owned part of the grid */
89: for (i = info->xs; i < info->xs+info->xm; ++i) {
90: if (i == 0) {
91: f[i].s = u[i].s - user->sl;
92: f[i].v = u[i].v - user->vl;
93: f[i].p = u[i].p - user->pl;
94: } else {
95: PetscScalar K = 2*dx/(dx/Kappa[i] + dx/Kappa[i-1]);
96: PetscReal lambdaWet = kappaWet*pow(u[i].s, alpha);
97: PetscReal lambda = lambdaWet + kappaNoWet*pow(1-u[i].s, beta);
98: PetscReal lambdaWetL = kappaWet*pow(u[i-1].s, alpha);
99: PetscReal lambdaL = lambdaWetL + kappaNoWet*pow(1-u[i-1].s, beta);
101: f[i].s = phi*(u[i].s - uold[i].s) + (dt/dx)*((lambdaWet/lambda)*u[i].v - (lambdaWetL/lambdaL)*u[i-1].v);
103: f[i].v = u[i].v + K*lambda*(u[i].p - u[i-1].p)/dx;
105: /*pxx = (2.0*u[i].p - u[i-1].p - u[i+1].p)/dx;*/
106: f[i].p = u[i].v - u[i-1].v;
107: }
108: }
109: DMDAVecRestoreArray(info->da, user->uold, &uold);
110: DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa);
111: /* PetscLogFlops(11.0*info->ym*info->xm); */
113: DMRestoreGlobalVector(user->cda, &L);
114: return(0);
115: }
119: int main(int argc, char **argv)
120: {
121: SNES snes; /* nonlinear solver */
122: DM da; /* grid */
123: Vec u; /* solution vector */
124: AppCtx user; /* user-defined work context */
125: PetscReal t; /* time */
127: PetscInt n;
129: PetscInitialize(&argc, &argv, NULL, help);
130: /* Create solver */
131: SNESCreate(PETSC_COMM_WORLD, &snes);
132: /* Create mesh */
133: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-4,3,1,NULL,&da);
134: DMSetApplicationContext(da, &user);
135: SNESSetDM(snes, da);
136: /* Create coefficient */
137: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-4,1,1,NULL,&user.cda);
138: DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
139: DMGetGlobalVector(user.cda, &user.Kappa);
140: FormPermeability(user.cda, user.Kappa, &user);
141: /* Setup Problem */
142: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
143: DMGetGlobalVector(da, &u);
144: DMGetGlobalVector(da, &user.uold);
146: user.sl = 1.0;
147: user.vl = 0.1;
148: user.pl = 1.0;
149: user.phi = 1.0;
151: user.kappaWet = 1.0;
152: user.kappaNoWet = 0.3;
154: /* Time Loop */
155: user.dt = 0.1;
156: for (n = 0; n < 100; ++n, t += user.dt) {
157: PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", t);
158: VecView(u, PETSC_VIEWER_DRAW_WORLD);
159: /* Solve */
160: SNESSetFromOptions(snes);
161: SNESSolve(snes, NULL, u);
162: /* Update */
163: VecCopy(u, user.uold);
165: VecView(u, PETSC_VIEWER_DRAW_WORLD);
166: }
167: /* Cleanup */
168: DMRestoreGlobalVector(da, &u);
169: DMRestoreGlobalVector(da, &user.uold);
170: DMRestoreGlobalVector(user.cda, &user.Kappa);
171: DMDestroy(&user.cda);
172: DMDestroy(&da);
173: SNESDestroy(&snes);
174: PetscFinalize();
175: return 0;
176: }