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