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(0); 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(0); 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 PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 122 /* Create solver */ 123 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 124 /* Create mesh */ 125 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,3,1,NULL,&da)); 126 PetscCall(DMSetFromOptions(da)); 127 PetscCall(DMSetUp(da)); 128 PetscCall(DMSetApplicationContext(da, &user)); 129 PetscCall(SNESSetDM(snes, da)); 130 /* Create coefficient */ 131 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,1,1,NULL,&user.cda)); 132 PetscCall(DMSetFromOptions(user.cda)); 133 PetscCall(DMSetUp(user.cda)); 134 PetscCall(DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 135 PetscCall(DMGetGlobalVector(user.cda, &user.Kappa)); 136 PetscCall(FormPermeability(user.cda, user.Kappa, &user)); 137 /* Setup Problem */ 138 PetscCall(DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user)); 139 PetscCall(DMGetGlobalVector(da, &u)); 140 PetscCall(DMGetGlobalVector(da, &user.uold)); 141 142 user.sl = 1.0; 143 user.vl = 0.1; 144 user.pl = 1.0; 145 user.phi = 1.0; 146 147 user.kappaWet = 1.0; 148 user.kappaNoWet = 0.3; 149 150 /* Time Loop */ 151 user.dt = 0.1; 152 for (n = 0; n < 100; ++n, t += user.dt) { 153 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", (double)t)); 154 PetscCall(VecView(u, PETSC_VIEWER_DRAW_WORLD)); 155 /* Solve */ 156 PetscCall(SNESSetFromOptions(snes)); 157 PetscCall(SNESSolve(snes, NULL, u)); 158 /* Update */ 159 PetscCall(VecCopy(u, user.uold)); 160 161 PetscCall(VecView(u, PETSC_VIEWER_DRAW_WORLD)); 162 } 163 /* Cleanup */ 164 PetscCall(DMRestoreGlobalVector(da, &u)); 165 PetscCall(DMRestoreGlobalVector(da, &user.uold)); 166 PetscCall(DMRestoreGlobalVector(user.cda, &user.Kappa)); 167 PetscCall(DMDestroy(&user.cda)); 168 PetscCall(DMDestroy(&da)); 169 PetscCall(SNESDestroy(&snes)); 170 PetscCall(PetscFinalize()); 171 return 0; 172 } 173 174 /*TEST 175 176 test: 177 suffix: 0 178 requires: !single 179 args: -snes_converged_reason -snes_monitor_short 180 181 TEST*/ 182