Actual source code: ex42.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: static char help[] = "Meinhard't activator-inhibitor model to test TS domain error feature.\n";

  3: /*
  4:    The activator-inhibitor on a line is described by the PDE:

  6:    da/dt = \alpha a^2 / (1 + \beta h) + \rho_a - \mu_a a + D_a d^2 a/ dx^2
  7:    dh/dt = \alpha a^2 + \rho_h - \mu_h h + D_h d^2 h/ dx^2

  9:    The PDE part will be solve by finite-difference on the line of cells.
 10:  */

 12: #include <petscts.h>

 14: typedef struct {
 15:   PetscInt  nb_cells;
 16:   PetscReal alpha;
 17:   PetscReal beta;
 18:   PetscReal rho_a;
 19:   PetscReal rho_h;
 20:   PetscReal mu_a;
 21:   PetscReal mu_h;
 22:   PetscReal D_a;
 23:   PetscReal D_h;
 24: } AppCtx;

 28: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec DXDT, void* ptr)
 29: {
 30:   AppCtx*           user = (AppCtx*)ptr;
 31:   PetscInt          nb_cells, i;
 32:   PetscReal         alpha, beta;
 33:   PetscReal         rho_a, mu_a, D_a;
 34:   PetscReal         rho_h, mu_h, D_h;
 35:   PetscReal         a, h, da, dh, d2a, d2h;
 36:   PetscErrorCode    ierr;
 37:   PetscScalar       *dxdt;
 38:   const PetscScalar *x;

 41:   nb_cells = user->nb_cells;
 42:   alpha    = user->alpha;
 43:   beta     = user->beta;
 44:   rho_a    = user->rho_a;
 45:   mu_a     = user->mu_a;
 46:   D_a      = user->D_a;
 47:   rho_h    = user->rho_h;
 48:   mu_h     = user->mu_h;
 49:   D_h      = user->D_h;

 51:   VecGetArrayRead(X, &x);
 52:   VecGetArray(DXDT, &dxdt);

 54:   for(i = 0 ; i < nb_cells ; i++) {
 55:     a = x[2*i];
 56:     h = x[2*i+1];
 57:     // Reaction:
 58:     da = alpha * a*a / (1. + beta * h) + rho_a - mu_a * a;
 59:     dh = alpha * a*a + rho_h - mu_h*h;
 60:     // Diffusion:
 61:     d2a = d2h = 0.;
 62:     if(i > 0) {
 63:       d2a += (x[2*(i-1)] - a);
 64:       d2h += (x[2*(i-1)+1] - h);
 65:     }
 66:     if(i < nb_cells-1) {
 67:       d2a += (x[2*(i+1)] - a);
 68:       d2h += (x[2*(i+1)+1] - h);
 69:     }
 70:     dxdt[2*i] = da + D_a*d2a;
 71:     dxdt[2*i+1] = dh + D_h*d2h;
 72:   }
 73:   VecRestoreArray(DXDT, &dxdt);
 74:   VecRestoreArrayRead(X, &x);
 75:   return(0);
 76: }

 80: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
 81: {
 82:   AppCtx            *user = (AppCtx*)ptr;
 83:   PetscInt          nb_cells, i, idx;
 84:   PetscReal         alpha, beta;
 85:   PetscReal         mu_a, D_a;
 86:   PetscReal         mu_h, D_h;
 87:   PetscReal         a, h;
 88:   const PetscScalar *x;
 89:   PetscScalar       va[4], vh[4];
 90:   PetscInt          ca[4], ch[4], rowa, rowh;
 91:   PetscErrorCode    ierr;

 94:   nb_cells = user->nb_cells;
 95:   alpha    = user->alpha;
 96:   beta     = user->beta;
 97:   mu_a     = user->mu_a;
 98:   D_a      = user->D_a;
 99:   mu_h     = user->mu_h;
100:   D_h      = user->D_h;

102:   VecGetArrayRead(X, &x);
103:   for(i = 0; i < nb_cells ; ++i) {
104:     rowa = 2*i;
105:     rowh = 2*i+1;
106:     a = x[2*i];
107:     h = x[2*i+1];
108:     ca[0] = ch[1] = 2*i;
109:     va[0] = 2*alpha*a / (1.+beta*h) - mu_a;
110:     vh[1] = 2*alpha*a;
111:     ca[1] = ch[0] = 2*i+1;
112:     va[1] = -beta*alpha*a*a / ((1.+beta*h)*(1.+beta*h));
113:     vh[0] = -mu_h;
114:     idx = 2;
115:     if(i > 0) {
116:       ca[idx] = 2*(i-1);
117:       ch[idx] = 2*(i-1)+1;
118:       va[idx] = D_a;
119:       vh[idx] = D_h;
120:       va[0] -= D_a;
121:       vh[0] -= D_h;
122:       idx++;
123:     }
124:     if(i < nb_cells-1) {
125:       ca[idx] = 2*(i+1);
126:       ch[idx] = 2*(i+1)+1;
127:       va[idx] = D_a;
128:       vh[idx] = D_h;
129:       va[0] -= D_a;
130:       vh[0] -= D_h;
131:       idx++;
132:     }
133:     MatSetValues(B, 1, &rowa, idx, ca, va, INSERT_VALUES);
134:     MatSetValues(B, 1, &rowh, idx, ch, vh, INSERT_VALUES);
135:   }
136:   VecRestoreArrayRead(X, &x);
137:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
139:   if (J != B) {
140:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
141:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
142:   }
143:   return(0);
144: }

148: PetscErrorCode DomainErrorFunction(TS ts, PetscReal t, Vec Y, PetscBool *accept)
149: {
150:   AppCtx            *user;
151:   PetscReal         dt;
152:   PetscErrorCode    ierr;
153:   const PetscScalar *x;
154:   PetscInt          nb_cells, i;

156:   TSGetApplicationContext(ts, &user);
157:   nb_cells = user->nb_cells;
158:   VecGetArrayRead(Y, &x);
159:   for(i = 0 ; i < 2*nb_cells ; ++i) {
160:     if(PetscRealPart(x[i]) < 0) {
161:       TSGetTimeStep(ts, &dt);
162:       PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t);
163:       *accept = PETSC_FALSE;
164:       break;
165:     }
166:   }
167:   VecRestoreArrayRead(Y, &x);
168:   return(0);
169: }

173: PetscErrorCode FormInitialState(Vec X, AppCtx* user)
174: {
176:   PetscRandom    R;

179:   PetscRandomCreate(PETSC_COMM_WORLD, &R);
180:   PetscRandomSetFromOptions(R);
181:   PetscRandomSetInterval(R, 0., 10.);

183:   /*
184:    * Initialize the state vector
185:    */
186:   VecSetRandom(X, R);
187:   PetscRandomDestroy(&R);
188:   return(0);
189: }

193: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
194: {
195:   PetscErrorCode    ierr;
196:   const PetscScalar *x;
197:   PetscInt          i;
198:   PetscInt          nb_cells = user->nb_cells;

201:   VecGetArrayRead(X, &x);
202:   PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");
203:   for(i = 0 ; i < nb_cells ; i++) {
204:     PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2*i], (double)x[2*i+1]);
205:   }
206:   VecRestoreArrayRead(X, &x);
207:   return(0);
208: }

212: int main(int argc, char **argv)
213: {
214:   TS             ts;       /* time-stepping context */
215:   Vec            x;       /* State vector */
216:   Mat            J; /* Jacobian matrix */
217:   AppCtx         user; /* user-defined context */
219:   PetscReal      ftime;
220:   PetscInt       its;
221:   PetscMPIInt    size;

223:   PetscInitialize(&argc, &argv, NULL, help);
224:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
225:   if(size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only");

227:   /*
228:    * Allow user to set the grid dimensions and the equations parameters
229:    */

231:   user.nb_cells = 50;
232:   user.alpha = 10.;
233:   user.beta = 1.;
234:   user.rho_a = 1.;
235:   user.rho_h = 2.;
236:   user.mu_a = 2.;
237:   user.mu_h = 3.;
238:   user.D_a = 0.;
239:   user.D_h = 30.;

241:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem settings", "PROBLEM");
242:   PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c",user.nb_cells, &user.nb_cells,NULL);
243:   PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c",user.alpha, &user.alpha,NULL);
244:   PetscOptionsReal("-beta", "Inhibition factor", "ex42.c",user.beta, &user.beta,NULL);
245:   PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c",user.rho_a, &user.rho_a,NULL);
246:   PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c",user.mu_a, &user.mu_a,NULL);
247:   PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c",user.D_a, &user.D_a,NULL);
248:   PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c",user.rho_h, &user.rho_h,NULL);
249:   PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c",user.mu_h, &user.mu_h,NULL);
250:   PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c",user.D_h, &user.D_h,NULL);
251:   PetscOptionsEnd();

253:   PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %D\n", user.nb_cells);
254:   PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", user.alpha);
255:   PetscPrintf(PETSC_COMM_WORLD, "beta:  %5.5g\n", user.beta);
256:   PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", user.rho_a);
257:   PetscPrintf(PETSC_COMM_WORLD, "mu_a:  %5.5g\n", user.mu_a);
258:   PetscPrintf(PETSC_COMM_WORLD, "D_a:   %5.5g\n", user.D_a);
259:   PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", user.rho_h);
260:   PetscPrintf(PETSC_COMM_WORLD, "mu_h:  %5.5g\n", user.mu_h);
261:   PetscPrintf(PETSC_COMM_WORLD, "D_h:   %5.5g\n", user.D_h);

263:   /*
264:    * Create vector to hold the solution
265:    */
266:   VecCreateSeq(PETSC_COMM_WORLD, 2*user.nb_cells, &x);

268:   /*
269:    * Create time-stepper context
270:    */
271:   TSCreate(PETSC_COMM_WORLD, &ts);
272:   TSSetProblemType(ts, TS_NONLINEAR);

274:   /*
275:    * Tell the time-stepper context where to compute the solution
276:    */
277:   TSSetSolution(ts, x);

279:   /*
280:    * Allocate the jacobian matrix
281:    */
282:   MatCreateSeqAIJ(PETSC_COMM_WORLD, 2*user.nb_cells, 2*user.nb_cells, 4, 0, &J);

284:   /*
285:    * Provide the call-back for the non-linear function we are evaluating.
286:    */
287:   TSSetRHSFunction(ts, NULL, RHSFunction, &user);

289:   /*
290:    * Set the Jacobian matrix and the function user to compute Jacobians
291:    */
292:   TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);

294:   /*
295:    * Set the function checking the domain
296:    */
297:   TSSetFunctionDomainError(ts, &DomainErrorFunction);

299:   /*
300:    * Initialize the problem with random values
301:    */
302:   FormInitialState(x, &user);

304:   /*
305:    * Read the solver type from options
306:    */
307:   TSSetType(ts, TSPSEUDO);

309:   /*
310:    * Set a large number of timesteps and final duration time to insure
311:    * convergenge to steady state
312:    */
313:   TSSetDuration(ts, 5000, 1e12);
314:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

316:   /*
317:    * Set a larger number of potential errors
318:    */
319:   TSSetMaxStepRejections(ts, 50);

321:   /*
322:    * Also start with a very small dt
323:    */
324:   TSSetTimeStep(ts, 0.05);

326:   /*
327:    * Set a larger time step increment
328:    */
329:   TSPseudoSetTimeStepIncrement(ts, 1.5);

331:   /*
332:    * Let the user personalise TS
333:    */
334:   TSSetFromOptions(ts);

336:   /*
337:    * Set the context for the time stepper
338:    */
339:   TSSetApplicationContext(ts, &user);

341:   /*
342:    * Setup the time stepper, ready for evaluation
343:    */
344:   TSSetUp(ts);

346:   /*
347:    * Perform the solve.
348:    */
349:   TSSolve(ts, x);
350:   TSGetSolveTime(ts, &ftime);
351:   TSGetTimeStepNumber(ts,&its);
352:   PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %D, final time: %4.2e\nResult:\n\n", its, (double)ftime);
353:   PrintSolution(x, &user);

355:   /*
356:    * Free the data structures
357:    */
358:   VecDestroy(&x);
359:   MatDestroy(&J);
360:   TSDestroy(&ts);
361:   PetscFinalize();
362:     return 0;
363: }