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 Set initial conditions 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 ierr = InitialConditions(da,U);CHKERRQ(ierr); 98 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Set solver options 102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 103 ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr); 104 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 105 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 106 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 107 108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 Solve nonlinear system 110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 111 ierr = TSSolve(ts,U);CHKERRQ(ierr); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Free work space. All PETSc objects should be destroyed when they 115 are no longer needed. 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 ierr = VecDestroy(&U);CHKERRQ(ierr); 118 ierr = TSDestroy(&ts);CHKERRQ(ierr); 119 ierr = DMDestroy(&da);CHKERRQ(ierr); 120 121 ierr = PetscFinalize(); 122 return ierr; 123 } 124 /* ------------------------------------------------------------------- */ 125 /* 126 IFunction - Evaluates nonlinear function, F(U). 127 128 Input Parameters: 129 . ts - the TS context 130 . U - input vector 131 . ptr - optional user-defined context, as set by SNESSetFunction() 132 133 Output Parameter: 134 . F - function vector 135 */ 136 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr) 137 { 138 AppCtx *appctx = (AppCtx*)ptr; 139 DM da; 140 PetscErrorCode ierr; 141 PetscInt i,Mx,xs,xm; 142 PetscReal hx,sx; 143 PetscScalar rho,c,rhoxx,cxx,cx,rhox,kcxrhox; 144 Field *u,*f,*udot; 145 Vec localU; 146 147 PetscFunctionBegin; 148 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 149 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 150 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); 151 152 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 153 154 /* 155 Scatter ghost points to local vector,using the 2-step process 156 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 157 By placing code between these two statements, computations can be 158 done while messages are in transition. 159 */ 160 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 161 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 162 163 /* 164 Get pointers to vector data 165 */ 166 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 167 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr); 168 ierr = DMDAVecGetArrayWrite(da,F,&f);CHKERRQ(ierr); 169 170 /* 171 Get local grid boundaries 172 */ 173 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 174 175 if (!xs) { 176 f[0].rho = udot[0].rho; /* u[0].rho - 0.0; */ 177 f[0].c = udot[0].c; /* u[0].c - 1.0; */ 178 xs++; 179 xm--; 180 } 181 if (xs+xm == Mx) { 182 f[Mx-1].rho = udot[Mx-1].rho; /* u[Mx-1].rho - 1.0; */ 183 f[Mx-1].c = udot[Mx-1].c; /* u[Mx-1].c - 0.0; */ 184 xm--; 185 } 186 187 /* 188 Compute function over the locally owned part of the grid 189 */ 190 for (i=xs; i<xs+xm; i++) { 191 rho = u[i].rho; 192 rhoxx = (-2.0*rho + u[i-1].rho + u[i+1].rho)*sx; 193 c = u[i].c; 194 cxx = (-2.0*c + u[i-1].c + u[i+1].c)*sx; 195 196 if (!appctx->upwind) { 197 rhox = .5*(u[i+1].rho - u[i-1].rho)/hx; 198 cx = .5*(u[i+1].c - u[i-1].c)/hx; 199 kcxrhox = appctx->kappa*(cxx*rho + cx*rhox); 200 } else { 201 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; 202 } 203 204 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; 205 f[i].c = udot[i].c - appctx->delta*cxx + appctx->lambda*c + appctx->alpha*rho*c/(appctx->gamma + c); 206 } 207 208 /* 209 Restore vectors 210 */ 211 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 212 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr); 213 ierr = DMDAVecRestoreArrayWrite(da,F,&f);CHKERRQ(ierr); 214 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 /* ------------------------------------------------------------------- */ 219 PetscErrorCode InitialConditions(DM da,Vec U) 220 { 221 PetscErrorCode ierr; 222 PetscInt i,xs,xm,Mx; 223 Field *u; 224 PetscReal hx,x; 225 226 PetscFunctionBegin; 227 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); 228 229 hx = 1.0/(PetscReal)(Mx-1); 230 231 /* 232 Get pointers to vector data 233 */ 234 ierr = DMDAVecGetArrayWrite(da,U,&u);CHKERRQ(ierr); 235 236 /* 237 Get local grid boundaries 238 */ 239 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 240 241 /* 242 Compute function over the locally owned part of the grid 243 */ 244 for (i=xs; i<xs+xm; i++) { 245 x = i*hx; 246 if (i < Mx-1) u[i].rho = 0.0; 247 else u[i].rho = 1.0; 248 u[i].c = PetscCosReal(.5*PETSC_PI*x); 249 } 250 251 /* 252 Restore vectors 253 */ 254 ierr = DMDAVecRestoreArrayWrite(da,U,&u);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 /*TEST 259 260 test: 261 args: -pc_type lu -da_refine 2 -ts_view -ts_monitor -ts_max_time 1 262 requires: double 263 264 test: 265 suffix: 2 266 args: -pc_type lu -da_refine 2 -ts_view -ts_monitor_draw_solution -ts_monitor -ts_max_time 1 267 requires: x double 268 output_file: output/ex4_1.out 269 270 TEST*/ 271