Actual source code: ex42.c

  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;

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

 38:   nb_cells = user->nb_cells;
 39:   alpha    = user->alpha;
 40:   beta     = user->beta;
 41:   rho_a    = user->rho_a;
 42:   mu_a     = user->mu_a;
 43:   D_a      = user->D_a;
 44:   rho_h    = user->rho_h;
 45:   mu_h     = user->mu_h;
 46:   D_h      = user->D_h;

 48:   VecGetArrayRead(X, &x);
 49:   VecGetArray(DXDT, &dxdt);

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

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

 88:   nb_cells = user->nb_cells;
 89:   alpha    = user->alpha;
 90:   beta     = user->beta;
 91:   mu_a     = user->mu_a;
 92:   D_a      = user->D_a;
 93:   mu_h     = user->mu_h;
 94:   D_h      = user->D_h;

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

140: PetscErrorCode DomainErrorFunction(TS ts, PetscReal t, Vec Y, PetscBool *accept)
141: {
142:   AppCtx            *user;
143:   PetscReal          dt;
144:   const PetscScalar *x;
145:   PetscInt           nb_cells, i;

148:   TSGetApplicationContext(ts, &user);
149:   nb_cells = user->nb_cells;
150:   VecGetArrayRead(Y, &x);
151:   for (i = 0; i < 2 * nb_cells; ++i) {
152:     if (PetscRealPart(x[i]) < 0) {
153:       TSGetTimeStep(ts, &dt);
154:       PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t);
155:       *accept = PETSC_FALSE;
156:       break;
157:     }
158:   }
159:   VecRestoreArrayRead(Y, &x);
160:   return 0;
161: }

163: PetscErrorCode FormInitialState(Vec X, AppCtx *user)
164: {
165:   PetscRandom R;

168:   PetscRandomCreate(PETSC_COMM_WORLD, &R);
169:   PetscRandomSetFromOptions(R);
170:   PetscRandomSetInterval(R, 0., 10.);

172:   /*
173:    * Initialize the state vector
174:    */
175:   VecSetRandom(X, R);
176:   PetscRandomDestroy(&R);
177:   return 0;
178: }

180: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
181: {
182:   const PetscScalar *x;
183:   PetscInt           i;
184:   PetscInt           nb_cells = user->nb_cells;

187:   VecGetArrayRead(X, &x);
188:   PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");
189:   for (i = 0; i < nb_cells; i++) PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2 * i], (double)x[2 * i + 1]);
190:   VecRestoreArrayRead(X, &x);
191:   return 0;
192: }

194: int main(int argc, char **argv)
195: {
196:   TS          ts;   /* time-stepping context */
197:   Vec         x;    /* State vector */
198:   Mat         J;    /* Jacobian matrix */
199:   AppCtx      user; /* user-defined context */
200:   PetscReal   ftime;
201:   PetscInt    its;
202:   PetscMPIInt size;

205:   PetscInitialize(&argc, &argv, NULL, help);
206:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

209:   /*
210:    * Allow user to set the grid dimensions and the equations parameters
211:    */

213:   user.nb_cells = 50;
214:   user.alpha    = 10.;
215:   user.beta     = 1.;
216:   user.rho_a    = 1.;
217:   user.rho_h    = 2.;
218:   user.mu_a     = 2.;
219:   user.mu_h     = 3.;
220:   user.D_a      = 0.;
221:   user.D_h      = 30.;

223:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem settings", "PROBLEM");
224:   PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c", user.nb_cells, &user.nb_cells, NULL);
225:   PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c", user.alpha, &user.alpha, NULL);
226:   PetscOptionsReal("-beta", "Inhibition factor", "ex42.c", user.beta, &user.beta, NULL);
227:   PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c", user.rho_a, &user.rho_a, NULL);
228:   PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c", user.mu_a, &user.mu_a, NULL);
229:   PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c", user.D_a, &user.D_a, NULL);
230:   PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c", user.rho_h, &user.rho_h, NULL);
231:   PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c", user.mu_h, &user.mu_h, NULL);
232:   PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c", user.D_h, &user.D_h, NULL);
233:   PetscOptionsEnd();

235:   PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %" PetscInt_FMT "\n", user.nb_cells);
236:   PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", (double)user.alpha);
237:   PetscPrintf(PETSC_COMM_WORLD, "beta:  %5.5g\n", (double)user.beta);
238:   PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", (double)user.rho_a);
239:   PetscPrintf(PETSC_COMM_WORLD, "mu_a:  %5.5g\n", (double)user.mu_a);
240:   PetscPrintf(PETSC_COMM_WORLD, "D_a:   %5.5g\n", (double)user.D_a);
241:   PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", (double)user.rho_h);
242:   PetscPrintf(PETSC_COMM_WORLD, "mu_h:  %5.5g\n", (double)user.mu_h);
243:   PetscPrintf(PETSC_COMM_WORLD, "D_h:   %5.5g\n", (double)user.D_h);

245:   /*
246:    * Create vector to hold the solution
247:    */
248:   VecCreateSeq(PETSC_COMM_WORLD, 2 * user.nb_cells, &x);

250:   /*
251:    * Create time-stepper context
252:    */
253:   TSCreate(PETSC_COMM_WORLD, &ts);
254:   TSSetProblemType(ts, TS_NONLINEAR);

256:   /*
257:    * Tell the time-stepper context where to compute the solution
258:    */
259:   TSSetSolution(ts, x);

261:   /*
262:    * Allocate the jacobian matrix
263:    */
264:   MatCreateSeqAIJ(PETSC_COMM_WORLD, 2 * user.nb_cells, 2 * user.nb_cells, 4, 0, &J);

266:   /*
267:    * Provide the call-back for the non-linear function we are evaluating.
268:    */
269:   TSSetRHSFunction(ts, NULL, RHSFunction, &user);

271:   /*
272:    * Set the Jacobian matrix and the function user to compute Jacobians
273:    */
274:   TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);

276:   /*
277:    * Set the function checking the domain
278:    */
279:   TSSetFunctionDomainError(ts, &DomainErrorFunction);

281:   /*
282:    * Initialize the problem with random values
283:    */
284:   FormInitialState(x, &user);

286:   /*
287:    * Read the solver type from options
288:    */
289:   TSSetType(ts, TSPSEUDO);

291:   /*
292:    * Set a large number of timesteps and final duration time to insure
293:    * convergenge to steady state
294:    */
295:   TSSetMaxSteps(ts, 2147483647);
296:   TSSetMaxTime(ts, 1.e12);
297:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

299:   /*
300:    * Set a larger number of potential errors
301:    */
302:   TSSetMaxStepRejections(ts, 50);

304:   /*
305:    * Also start with a very small dt
306:    */
307:   TSSetTimeStep(ts, 0.05);

309:   /*
310:    * Set a larger time step increment
311:    */
312:   TSPseudoSetTimeStepIncrement(ts, 1.5);

314:   /*
315:    * Let the user personalise TS
316:    */
317:   TSSetFromOptions(ts);

319:   /*
320:    * Set the context for the time stepper
321:    */
322:   TSSetApplicationContext(ts, &user);

324:   /*
325:    * Setup the time stepper, ready for evaluation
326:    */
327:   TSSetUp(ts);

329:   /*
330:    * Perform the solve.
331:    */
332:   TSSolve(ts, x);
333:   TSGetSolveTime(ts, &ftime);
334:   TSGetStepNumber(ts, &its);
335:   PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %" PetscInt_FMT ", final time: %4.2e\nResult:\n\n", its, (double)ftime);
336:   PrintSolution(x, &user);

338:   /*
339:    * Free the data structures
340:    */
341:   VecDestroy(&x);
342:   MatDestroy(&J);
343:   TSDestroy(&ts);
344:   PetscFinalize();
345:   return 0;
346: }

348: /*TEST
349:     build:
350:       requires: !single !complex

352:     test:
353:       args: -ts_max_steps 8
354:       output_file: output/ex42.out

356: TEST*/