1 2 static char help[] = "Chemo-taxis Problems from Mathematical Biology.\n"; 3 4 /* 5 Page 18, Chemo-taxis Problems from Mathematical Biology 6 7 rho_t = 8 c_t = 9 10 Further discusson on Page 134 and in 2d on Page 409 11 */ 12 13 /* 14 15 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 16 Include "petscts.h" so that we can use SNES solvers. Note that this 17 file automatically includes: 18 petscsys.h - base PETSc routines petscvec.h - vectors 19 petscmat.h - matrices 20 petscis.h - index sets petscksp.h - Krylov subspace methods 21 petscviewer.h - viewers petscpc.h - preconditioners 22 petscksp.h - linear solvers 23 */ 24 25 #include <petscdm.h> 26 #include <petscdmda.h> 27 #include <petscts.h> 28 29 typedef struct { 30 PetscScalar rho,c; 31 } Field; 32 33 typedef struct { 34 PetscScalar epsilon,delta,alpha,beta,gamma,kappa,lambda,mu,cstar; 35 PetscBool upwind; 36 } AppCtx; 37 38 /* 39 User-defined routines 40 */ 41 extern PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*),InitialConditions(DM,Vec); 42 43 int main(int argc,char **argv) 44 { 45 TS ts; /* nonlinear solver */ 46 Vec U; /* solution, residual vectors */ 47 PetscErrorCode ierr; 48 DM da; 49 AppCtx appctx; 50 51 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 52 Initialize program 53 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 54 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 55 56 appctx.epsilon = 1.0e-3; 57 appctx.delta = 1.0; 58 appctx.alpha = 10.0; 59 appctx.beta = 4.0; 60 appctx.gamma = 1.0; 61 appctx.kappa = .75; 62 appctx.lambda = 1.0; 63 appctx.mu = 100.; 64 appctx.cstar = .2; 65 appctx.upwind = PETSC_TRUE; 66 67 ierr = PetscOptionsGetScalar(NULL,NULL,"-delta",&appctx.delta,NULL);CHKERRQ(ierr); 68 ierr = PetscOptionsGetBool(NULL,NULL,"-upwind",&appctx.upwind,NULL);CHKERRQ(ierr); 69 70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 71 Create distributed array (DMDA) to manage parallel grid and vectors 72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 73 ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,8,2,1,NULL,&da);CHKERRQ(ierr); 74 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 75 ierr = DMSetUp(da);CHKERRQ(ierr); 76 ierr = DMDASetFieldName(da,0,"rho");CHKERRQ(ierr); 77 ierr = DMDASetFieldName(da,1,"c");CHKERRQ(ierr); 78 79 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 80 Extract global vectors from DMDA; then duplicate for remaining 81 vectors that are the same types 82 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 83 ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr); 84 85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86 Create timestepping solver context 87 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 88 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 89 ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 90 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 91 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 92 ierr = TSSetIFunction(ts,NULL,IFunction,&appctx);CHKERRQ(ierr); 93 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Set initial conditions 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 ierr = InitialConditions(da,U);CHKERRQ(ierr); 99 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Set solver options 103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104 ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr); 105 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 106 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 107 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Solve nonlinear system 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 ierr = TSSolve(ts,U);CHKERRQ(ierr); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Free work space. All PETSc objects should be destroyed when they 116 are no longer needed. 117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 118 ierr = VecDestroy(&U);CHKERRQ(ierr); 119 ierr = TSDestroy(&ts);CHKERRQ(ierr); 120 ierr = DMDestroy(&da);CHKERRQ(ierr); 121 122 ierr = PetscFinalize(); 123 return ierr; 124 } 125 /* ------------------------------------------------------------------- */ 126 /* 127 IFunction - Evaluates nonlinear function, F(U). 128 129 Input Parameters: 130 . ts - the TS context 131 . U - input vector 132 . ptr - optional user-defined context, as set by SNESSetFunction() 133 134 Output Parameter: 135 . F - function vector 136 */ 137 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr) 138 { 139 AppCtx *appctx = (AppCtx*)ptr; 140 DM da; 141 PetscErrorCode ierr; 142 PetscInt i,Mx,xs,xm; 143 PetscReal hx,sx; 144 PetscScalar rho,c,rhoxx,cxx,cx,rhox,kcxrhox; 145 Field *u,*f,*udot; 146 Vec localU; 147 148 PetscFunctionBegin; 149 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 150 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 151 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 152 153 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 154 155 /* 156 Scatter ghost points to local vector,using the 2-step process 157 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 158 By placing code between these two statements, computations can be 159 done while messages are in transition. 160 */ 161 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 162 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 163 164 /* 165 Get pointers to vector data 166 */ 167 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 168 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr); 169 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 170 171 /* 172 Get local grid boundaries 173 */ 174 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 175 176 if (!xs) { 177 f[0].rho = udot[0].rho; /* u[0].rho - 0.0; */ 178 f[0].c = udot[0].c; /* u[0].c - 1.0; */ 179 xs++; 180 xm--; 181 } 182 if (xs+xm == Mx) { 183 f[Mx-1].rho = udot[Mx-1].rho; /* u[Mx-1].rho - 1.0; */ 184 f[Mx-1].c = udot[Mx-1].c; /* u[Mx-1].c - 0.0; */ 185 xm--; 186 } 187 188 /* 189 Compute function over the locally owned part of the grid 190 */ 191 for (i=xs; i<xs+xm; i++) { 192 rho = u[i].rho; 193 rhoxx = (-2.0*rho + u[i-1].rho + u[i+1].rho)*sx; 194 c = u[i].c; 195 cxx = (-2.0*c + u[i-1].c + u[i+1].c)*sx; 196 197 if (!appctx->upwind) { 198 rhox = .5*(u[i+1].rho - u[i-1].rho)/hx; 199 cx = .5*(u[i+1].c - u[i-1].c)/hx; 200 kcxrhox = appctx->kappa*(cxx*rho + cx*rhox); 201 } else { 202 kcxrhox = appctx->kappa*((u[i+1].c - u[i].c)*u[i+1].rho - (u[i].c - u[i-1].c)*u[i].rho)*sx; 203 } 204 205 f[i].rho = udot[i].rho - appctx->epsilon*rhoxx + kcxrhox - appctx->mu*PetscAbsScalar(rho)*(1.0 - rho)*PetscMax(0,PetscRealPart(c - appctx->cstar)) + appctx->beta*rho; 206 f[i].c = udot[i].c - appctx->delta*cxx + appctx->lambda*c + appctx->alpha*rho*c/(appctx->gamma + c); 207 } 208 209 /* 210 Restore vectors 211 */ 212 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 213 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr); 214 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 215 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 /* ------------------------------------------------------------------- */ 220 PetscErrorCode InitialConditions(DM da,Vec U) 221 { 222 PetscErrorCode ierr; 223 PetscInt i,xs,xm,Mx; 224 Field *u; 225 PetscReal hx,x; 226 227 PetscFunctionBegin; 228 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 229 230 hx = 1.0/(PetscReal)(Mx-1); 231 232 /* 233 Get pointers to vector data 234 */ 235 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 236 237 /* 238 Get local grid boundaries 239 */ 240 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 241 242 /* 243 Compute function over the locally owned part of the grid 244 */ 245 for (i=xs; i<xs+xm; i++) { 246 x = i*hx; 247 if (x < 1.0) u[i].rho = 0.0; 248 else u[i].rho = 1.0; 249 u[i].c = PetscCosReal(.5*PETSC_PI*x); 250 } 251 252 /* 253 Restore vectors 254 */ 255 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 256 PetscFunctionReturn(0); 257 } 258 259 260 261 262 /*TEST 263 264 test: 265 args: -pc_type lu -da_refine 2 -ts_view -ts_monitor -ts_max_time 1 266 requires: double 267 268 test: 269 suffix: 2 270 args: -pc_type lu -da_refine 2 -ts_view -ts_monitor_draw_solution -ts_monitor -ts_max_time 1 271 requires: x double 272 output_file: output/ex4_1.out 273 274 TEST*/ 275