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 PetscErrorCode ierr; 41 42 PetscFunctionBeginUser; 43 ierr = DMGetCoordinateDM(da, &cda);CHKERRQ(ierr); 44 ierr = DMGetCoordinates(da, &c);CHKERRQ(ierr); 45 ierr = DMDAGetCorners(da, &xs,NULL,NULL, &xm,NULL,NULL);CHKERRQ(ierr); 46 ierr = DMDAVecGetArray(da, Kappa, &K);CHKERRQ(ierr); 47 ierr = DMDAVecGetArray(cda, c, &coords);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(da, Kappa, &K);CHKERRQ(ierr); 58 ierr = DMDAVecRestoreArray(cda, c, &coords);CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 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; 78 PetscErrorCode ierr; 79 80 PetscFunctionBeginUser; 81 ierr = DMGetGlobalVector(user->cda, &L);CHKERRQ(ierr); 82 83 ierr = DMDAVecGetArray(info->da, user->uold, &uold);CHKERRQ(ierr); 84 ierr = DMDAVecGetArray(user->cda, user->Kappa, &Kappa);CHKERRQ(ierr); 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)); 97 98 f[i].s = phi*(u[i].s - uold[i].s) + (dt/dx)*((lambdaWet/lambda)*u[i].v - (lambdaWetL/lambdaL)*u[i-1].v); 99 100 f[i].v = u[i].v + K*lambda*(u[i].p - u[i-1].p)/dx; 101 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 ierr = DMDAVecRestoreArray(info->da, user->uold, &uold);CHKERRQ(ierr); 107 ierr = DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa);CHKERRQ(ierr); 108 /* ierr = PetscLogFlops(11.0*info->ym*info->xm);CHKERRQ(ierr); */ 109 110 ierr = DMRestoreGlobalVector(user->cda, &L);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 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 */ 121 PetscErrorCode ierr; 122 PetscInt n; 123 124 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 125 /* Create solver */ 126 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 127 /* Create mesh */ 128 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,3,1,NULL,&da);CHKERRQ(ierr); 129 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 130 ierr = DMSetUp(da);CHKERRQ(ierr); 131 ierr = DMSetApplicationContext(da, &user);CHKERRQ(ierr); 132 ierr = SNESSetDM(snes, da);CHKERRQ(ierr); 133 /* Create coefficient */ 134 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,1,1,NULL,&user.cda);CHKERRQ(ierr); 135 ierr = DMSetFromOptions(user.cda);CHKERRQ(ierr); 136 ierr = DMSetUp(user.cda);CHKERRQ(ierr); 137 ierr = DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr); 138 ierr = DMGetGlobalVector(user.cda, &user.Kappa);CHKERRQ(ierr); 139 ierr = FormPermeability(user.cda, user.Kappa, &user);CHKERRQ(ierr); 140 /* Setup Problem */ 141 ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);CHKERRQ(ierr); 142 ierr = DMGetGlobalVector(da, &u);CHKERRQ(ierr); 143 ierr = DMGetGlobalVector(da, &user.uold);CHKERRQ(ierr); 144 145 user.sl = 1.0; 146 user.vl = 0.1; 147 user.pl = 1.0; 148 user.phi = 1.0; 149 150 user.kappaWet = 1.0; 151 user.kappaNoWet = 0.3; 152 153 /* Time Loop */ 154 user.dt = 0.1; 155 for (n = 0; n < 100; ++n, t += user.dt) { 156 ierr = PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", (double)t);CHKERRQ(ierr); 157 ierr = VecView(u, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 158 /* Solve */ 159 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 160 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 161 /* Update */ 162 ierr = VecCopy(u, user.uold);CHKERRQ(ierr); 163 164 ierr = VecView(u, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 165 } 166 /* Cleanup */ 167 ierr = DMRestoreGlobalVector(da, &u);CHKERRQ(ierr); 168 ierr = DMRestoreGlobalVector(da, &user.uold);CHKERRQ(ierr); 169 ierr = DMRestoreGlobalVector(user.cda, &user.Kappa);CHKERRQ(ierr); 170 ierr = DMDestroy(&user.cda);CHKERRQ(ierr); 171 ierr = DMDestroy(&da);CHKERRQ(ierr); 172 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 173 ierr = PetscFinalize(); 174 return ierr; 175 } 176 177 /*TEST 178 179 test: 180 suffix: 0 181 requires: !single 182 args: -snes_converged_reason -snes_monitor_short 183 184 TEST*/ 185