Actual source code: ex42.c

petsc-3.14.6 2021-03-30
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;

 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:   PetscErrorCode    ierr;
 35:   PetscScalar       *dxdt;
 36:   const PetscScalar *x;

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

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

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

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

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

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

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

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

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

171:   PetscRandomCreate(PETSC_COMM_WORLD, &R);
172:   PetscRandomSetFromOptions(R);
173:   PetscRandomSetInterval(R, 0., 10.);

175:   /*
176:    * Initialize the state vector
177:    */
178:   VecSetRandom(X, R);
179:   PetscRandomDestroy(&R);
180:   return(0);
181: }

183: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
184: {
185:   PetscErrorCode    ierr;
186:   const PetscScalar *x;
187:   PetscInt          i;
188:   PetscInt          nb_cells = user->nb_cells;

191:   VecGetArrayRead(X, &x);
192:   PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");
193:   for (i = 0 ; i < nb_cells ; i++) {
194:     PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2*i], (double)x[2*i+1]);
195:   }
196:   VecRestoreArrayRead(X, &x);
197:   return(0);
198: }

200: int main(int argc, char **argv)
201: {
202:   TS             ts;       /* time-stepping context */
203:   Vec            x;       /* State vector */
204:   Mat            J; /* Jacobian matrix */
205:   AppCtx         user; /* user-defined context */
207:   PetscReal      ftime;
208:   PetscInt       its;
209:   PetscMPIInt    size;

211:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
212:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
213:   if (size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only");

215:   /*
216:    * Allow user to set the grid dimensions and the equations parameters
217:    */

219:   user.nb_cells = 50;
220:   user.alpha = 10.;
221:   user.beta = 1.;
222:   user.rho_a = 1.;
223:   user.rho_h = 2.;
224:   user.mu_a = 2.;
225:   user.mu_h = 3.;
226:   user.D_a = 0.;
227:   user.D_h = 30.;

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

241:   PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %D\n", user.nb_cells);
242:   PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", (double)user.alpha);
243:   PetscPrintf(PETSC_COMM_WORLD, "beta:  %5.5g\n", (double)user.beta);
244:   PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", (double)user.rho_a);
245:   PetscPrintf(PETSC_COMM_WORLD, "mu_a:  %5.5g\n", (double)user.mu_a);
246:   PetscPrintf(PETSC_COMM_WORLD, "D_a:   %5.5g\n", (double)user.D_a);
247:   PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", (double)user.rho_h);
248:   PetscPrintf(PETSC_COMM_WORLD, "mu_h:  %5.5g\n", (double)user.mu_h);
249:   PetscPrintf(PETSC_COMM_WORLD, "D_h:   %5.5g\n", (double)user.D_h);

251:   /*
252:    * Create vector to hold the solution
253:    */
254:   VecCreateSeq(PETSC_COMM_WORLD, 2*user.nb_cells, &x);

256:   /*
257:    * Create time-stepper context
258:    */
259:   TSCreate(PETSC_COMM_WORLD, &ts);
260:   TSSetProblemType(ts, TS_NONLINEAR);

262:   /*
263:    * Tell the time-stepper context where to compute the solution
264:    */
265:   TSSetSolution(ts, x);

267:   /*
268:    * Allocate the jacobian matrix
269:    */
270:   MatCreateSeqAIJ(PETSC_COMM_WORLD, 2*user.nb_cells, 2*user.nb_cells, 4, 0, &J);

272:   /*
273:    * Provide the call-back for the non-linear function we are evaluating.
274:    */
275:   TSSetRHSFunction(ts, NULL, RHSFunction, &user);

277:   /*
278:    * Set the Jacobian matrix and the function user to compute Jacobians
279:    */
280:   TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);

282:   /*
283:    * Set the function checking the domain
284:    */
285:   TSSetFunctionDomainError(ts, &DomainErrorFunction);

287:   /*
288:    * Initialize the problem with random values
289:    */
290:   FormInitialState(x, &user);

292:   /*
293:    * Read the solver type from options
294:    */
295:   TSSetType(ts, TSPSEUDO);

297:   /*
298:    * Set a large number of timesteps and final duration time to insure
299:    * convergenge to steady state
300:    */
301:   TSSetMaxSteps(ts, 2147483647);
302:   TSSetMaxTime(ts, 1.e12);
303:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

305:   /*
306:    * Set a larger number of potential errors
307:    */
308:   TSSetMaxStepRejections(ts, 50);

310:   /*
311:    * Also start with a very small dt
312:    */
313:   TSSetTimeStep(ts, 0.05);

315:   /*
316:    * Set a larger time step increment
317:    */
318:   TSPseudoSetTimeStepIncrement(ts, 1.5);

320:   /*
321:    * Let the user personalise TS
322:    */
323:   TSSetFromOptions(ts);

325:   /*
326:    * Set the context for the time stepper
327:    */
328:   TSSetApplicationContext(ts, &user);

330:   /*
331:    * Setup the time stepper, ready for evaluation
332:    */
333:   TSSetUp(ts);

335:   /*
336:    * Perform the solve.
337:    */
338:   TSSolve(ts, x);
339:   TSGetSolveTime(ts, &ftime);
340:   TSGetStepNumber(ts,&its);
341:   PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %D, final time: %4.2e\nResult:\n\n", its, (double)ftime);
342:   PrintSolution(x, &user);

344:   /*
345:    * Free the data structures
346:    */
347:   VecDestroy(&x);
348:   MatDestroy(&J);
349:   TSDestroy(&ts);
350:   PetscFinalize();
351:   return ierr;
352: }

354: /*TEST
355:     build:
356:       requires: !single !complex

358:     test:
359:       args: -ts_max_steps 8
360:       output_file: output/ex42.out

362: TEST*/