Actual source code: ex4.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Chemo-taxis Problems from Mathematical Biology.\n";
4: /*
5: Page 18, Chemo-taxis Problems from Mathematical Biology
7: rho_t =
8: c_t =
10: Further discusson on Page 134 and in 2d on Page 409
11: */
13: /*
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: #include <petscdmda.h>
25: #include <petscts.h>
27: typedef struct {
28: PetscScalar rho,c;
29: } Field;
31: typedef struct {
32: PetscScalar epsilon,delta,alpha,beta,gamma,kappa,lambda,mu,cstar;
33: PetscBool upwind;
34: } AppCtx;
36: /*
37: User-defined routines
38: */
39: extern PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*),InitialConditions(DM,Vec);
43: int main(int argc,char **argv)
44: {
45: TS ts; /* nonlinear solver */
46: Vec U; /* solution, residual vectors */
48: DM da;
49: AppCtx appctx;
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Initialize program
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: PetscInitialize(&argc,&argv,(char*)0,help);
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;
67: PetscOptionsGetScalar(NULL,"-delta",&appctx.delta,NULL);
68: PetscOptionsGetBool(NULL,"-upwind",&appctx.upwind,NULL);
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Create distributed array (DMDA) to manage parallel grid and vectors
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: DMDACreate1d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE,-8,2,1,NULL,&da);
74: DMDASetFieldName(da,0,"rho");
75: DMDASetFieldName(da,1,"c");
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Extract global vectors from DMDA; then duplicate for remaining
79: vectors that are the same types
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: DMCreateGlobalVector(da,&U);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create timestepping solver context
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: TSCreate(PETSC_COMM_WORLD,&ts);
87: TSSetType(ts,TSROSW);
88: TSSetDM(ts,da);
89: TSSetProblemType(ts,TS_NONLINEAR);
90: TSSetIFunction(ts,NULL,IFunction,&appctx);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Set initial conditions
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: InitialConditions(da,U);
97: TSSetSolution(ts,U);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Set solver options
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: TSSetInitialTimeStep(ts,0.0,.0001);
103: TSSetDuration(ts,PETSC_DEFAULT,1.0);
104: TSSetFromOptions(ts);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Solve nonlinear system
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: TSSolve(ts,U);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Free work space. All PETSc objects should be destroyed when they
113: are no longer needed.
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: VecDestroy(&U);
116: TSDestroy(&ts);
117: DMDestroy(&da);
119: PetscFinalize();
120: return(0);
121: }
122: /* ------------------------------------------------------------------- */
125: /*
126: IFunction - Evaluates nonlinear function, F(U).
128: Input Parameters:
129: . ts - the TS context
130: . U - input vector
131: . ptr - optional user-defined context, as set by SNESSetFunction()
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;
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;
148: TSGetDM(ts,&da);
149: DMGetLocalVector(da,&localU);
150: 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);
152: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
161: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
163: /*
164: Get pointers to vector data
165: */
166: DMDAVecGetArray(da,localU,&u);
167: DMDAVecGetArray(da,Udot,&udot);
168: DMDAVecGetArray(da,F,&f);
170: /*
171: Get local grid boundaries
172: */
173: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
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: }
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;
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: }
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: }
208: /*
209: Restore vectors
210: */
211: DMDAVecRestoreArray(da,localU,&u);
212: DMDAVecRestoreArray(da,Udot,&udot);
213: DMDAVecRestoreArray(da,F,&f);
214: DMRestoreLocalVector(da,&localU);
215: return(0);
216: }
218: /* ------------------------------------------------------------------- */
221: PetscErrorCode InitialConditions(DM da,Vec U)
222: {
224: PetscInt i,xs,xm,Mx;
225: Field *u;
226: PetscReal hx,x;
229: 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);
231: hx = 1.0/(PetscReal)(Mx-1);
233: /*
234: Get pointers to vector data
235: */
236: DMDAVecGetArray(da,U,&u);
238: /*
239: Get local grid boundaries
240: */
241: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
243: /*
244: Compute function over the locally owned part of the grid
245: */
246: for (i=xs; i<xs+xm; i++) {
247: x = i*hx;
248: if (x < 1.0) u[i].rho = 0.0;
249: else u[i].rho = 1.0;
250: u[i].c = PetscCosScalar(.5*PETSC_PI*x);
251: }
253: /*
254: Restore vectors
255: */
256: DMDAVecRestoreArray(da,U,&u);
257: return(0);
258: }