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 { 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 */ 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 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