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;
37: nb_cells = user->nb_cells;
38: alpha = user->alpha;
39: beta = user->beta;
40: rho_a = user->rho_a;
41: mu_a = user->mu_a;
42: D_a = user->D_a;
43: rho_h = user->rho_h;
44: mu_h = user->mu_h;
45: D_h = user->D_h;
47: VecGetArrayRead(X, &x);
48: VecGetArray(DXDT, &dxdt);
50: for (i = 0 ; i < nb_cells ; i++) {
51: a = x[2*i];
52: h = x[2*i+1];
53: // Reaction:
54: da = alpha * a*a / (1. + beta * h) + rho_a - mu_a * a;
55: dh = alpha * a*a + rho_h - mu_h*h;
56: // Diffusion:
57: d2a = d2h = 0.;
58: if (i > 0) {
59: d2a += (x[2*(i-1)] - a);
60: d2h += (x[2*(i-1)+1] - h);
61: }
62: if (i < nb_cells-1) {
63: d2a += (x[2*(i+1)] - a);
64: d2h += (x[2*(i+1)+1] - h);
65: }
66: dxdt[2*i] = da + D_a*d2a;
67: dxdt[2*i+1] = dh + D_h*d2h;
68: }
69: VecRestoreArray(DXDT, &dxdt);
70: VecRestoreArrayRead(X, &x);
71: return 0;
72: }
74: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
75: {
76: AppCtx *user = (AppCtx*)ptr;
77: PetscInt nb_cells, i, idx;
78: PetscReal alpha, beta;
79: PetscReal mu_a, D_a;
80: PetscReal mu_h, D_h;
81: PetscReal a, h;
82: const PetscScalar *x;
83: PetscScalar va[4], vh[4];
84: PetscInt ca[4], ch[4], rowa, rowh;
86: nb_cells = user->nb_cells;
87: alpha = user->alpha;
88: beta = user->beta;
89: mu_a = user->mu_a;
90: D_a = user->D_a;
91: mu_h = user->mu_h;
92: D_h = user->D_h;
94: VecGetArrayRead(X, &x);
95: for (i = 0; i < nb_cells ; ++i) {
96: rowa = 2*i;
97: rowh = 2*i+1;
98: a = x[2*i];
99: h = x[2*i+1];
100: ca[0] = ch[1] = 2*i;
101: va[0] = 2*alpha*a / (1.+beta*h) - mu_a;
102: vh[1] = 2*alpha*a;
103: ca[1] = ch[0] = 2*i+1;
104: va[1] = -beta*alpha*a*a / ((1.+beta*h)*(1.+beta*h));
105: vh[0] = -mu_h;
106: idx = 2;
107: if (i > 0) {
108: ca[idx] = 2*(i-1);
109: ch[idx] = 2*(i-1)+1;
110: va[idx] = D_a;
111: vh[idx] = D_h;
112: va[0] -= D_a;
113: vh[0] -= D_h;
114: idx++;
115: }
116: if (i < nb_cells-1) {
117: ca[idx] = 2*(i+1);
118: ch[idx] = 2*(i+1)+1;
119: va[idx] = D_a;
120: vh[idx] = D_h;
121: va[0] -= D_a;
122: vh[0] -= D_h;
123: idx++;
124: }
125: MatSetValues(B, 1, &rowa, idx, ca, va, INSERT_VALUES);
126: MatSetValues(B, 1, &rowh, idx, ch, vh, INSERT_VALUES);
127: }
128: VecRestoreArrayRead(X, &x);
129: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
131: if (J != B) {
132: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
134: }
135: return 0;
136: }
138: PetscErrorCode DomainErrorFunction(TS ts, PetscReal t, Vec Y, PetscBool *accept)
139: {
140: AppCtx *user;
141: PetscReal dt;
142: const PetscScalar *x;
143: PetscInt nb_cells, i;
145: TSGetApplicationContext(ts, &user);
146: nb_cells = user->nb_cells;
147: VecGetArrayRead(Y, &x);
148: for (i = 0 ; i < 2*nb_cells ; ++i) {
149: if (PetscRealPart(x[i]) < 0) {
150: TSGetTimeStep(ts, &dt);
151: PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t);
152: *accept = PETSC_FALSE;
153: break;
154: }
155: }
156: VecRestoreArrayRead(Y, &x);
157: return 0;
158: }
160: PetscErrorCode FormInitialState(Vec X, AppCtx* user)
161: {
162: PetscRandom R;
164: PetscRandomCreate(PETSC_COMM_WORLD, &R);
165: PetscRandomSetFromOptions(R);
166: PetscRandomSetInterval(R, 0., 10.);
168: /*
169: * Initialize the state vector
170: */
171: VecSetRandom(X, R);
172: PetscRandomDestroy(&R);
173: return 0;
174: }
176: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
177: {
178: const PetscScalar *x;
179: PetscInt i;
180: PetscInt nb_cells = user->nb_cells;
182: VecGetArrayRead(X, &x);
183: PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");
184: for (i = 0 ; i < nb_cells ; i++) {
185: PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2*i], (double)x[2*i+1]);
186: }
187: VecRestoreArrayRead(X, &x);
188: return 0;
189: }
191: int main(int argc, char **argv)
192: {
193: TS ts; /* time-stepping context */
194: Vec x; /* State vector */
195: Mat J; /* Jacobian matrix */
196: AppCtx user; /* user-defined context */
198: PetscReal ftime;
199: PetscInt its;
200: PetscMPIInt size;
202: PetscInitialize(&argc, &argv, NULL, help);
203: MPI_Comm_size(PETSC_COMM_WORLD, &size);
206: /*
207: * Allow user to set the grid dimensions and the equations parameters
208: */
210: user.nb_cells = 50;
211: user.alpha = 10.;
212: user.beta = 1.;
213: user.rho_a = 1.;
214: user.rho_h = 2.;
215: user.mu_a = 2.;
216: user.mu_h = 3.;
217: user.D_a = 0.;
218: user.D_h = 30.;
220: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem settings", "PROBLEM");
221: PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c",user.nb_cells, &user.nb_cells,NULL);
222: PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c",user.alpha, &user.alpha,NULL);
223: PetscOptionsReal("-beta", "Inhibition factor", "ex42.c",user.beta, &user.beta,NULL);
224: PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c",user.rho_a, &user.rho_a,NULL);
225: PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c",user.mu_a, &user.mu_a,NULL);
226: PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c",user.D_a, &user.D_a,NULL);
227: PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c",user.rho_h, &user.rho_h,NULL);
228: PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c",user.mu_h, &user.mu_h,NULL);
229: PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c",user.D_h, &user.D_h,NULL);
230: PetscOptionsEnd();
232: PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %D\n", user.nb_cells);
233: PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", (double)user.alpha);
234: PetscPrintf(PETSC_COMM_WORLD, "beta: %5.5g\n", (double)user.beta);
235: PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", (double)user.rho_a);
236: PetscPrintf(PETSC_COMM_WORLD, "mu_a: %5.5g\n", (double)user.mu_a);
237: PetscPrintf(PETSC_COMM_WORLD, "D_a: %5.5g\n", (double)user.D_a);
238: PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", (double)user.rho_h);
239: PetscPrintf(PETSC_COMM_WORLD, "mu_h: %5.5g\n", (double)user.mu_h);
240: PetscPrintf(PETSC_COMM_WORLD, "D_h: %5.5g\n", (double)user.D_h);
242: /*
243: * Create vector to hold the solution
244: */
245: VecCreateSeq(PETSC_COMM_WORLD, 2*user.nb_cells, &x);
247: /*
248: * Create time-stepper context
249: */
250: TSCreate(PETSC_COMM_WORLD, &ts);
251: TSSetProblemType(ts, TS_NONLINEAR);
253: /*
254: * Tell the time-stepper context where to compute the solution
255: */
256: TSSetSolution(ts, x);
258: /*
259: * Allocate the jacobian matrix
260: */
261: MatCreateSeqAIJ(PETSC_COMM_WORLD, 2*user.nb_cells, 2*user.nb_cells, 4, 0, &J);
263: /*
264: * Provide the call-back for the non-linear function we are evaluating.
265: */
266: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
268: /*
269: * Set the Jacobian matrix and the function user to compute Jacobians
270: */
271: TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);
273: /*
274: * Set the function checking the domain
275: */
276: TSSetFunctionDomainError(ts, &DomainErrorFunction);
278: /*
279: * Initialize the problem with random values
280: */
281: FormInitialState(x, &user);
283: /*
284: * Read the solver type from options
285: */
286: TSSetType(ts, TSPSEUDO);
288: /*
289: * Set a large number of timesteps and final duration time to insure
290: * convergenge to steady state
291: */
292: TSSetMaxSteps(ts, 2147483647);
293: TSSetMaxTime(ts, 1.e12);
294: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
296: /*
297: * Set a larger number of potential errors
298: */
299: TSSetMaxStepRejections(ts, 50);
301: /*
302: * Also start with a very small dt
303: */
304: TSSetTimeStep(ts, 0.05);
306: /*
307: * Set a larger time step increment
308: */
309: TSPseudoSetTimeStepIncrement(ts, 1.5);
311: /*
312: * Let the user personalise TS
313: */
314: TSSetFromOptions(ts);
316: /*
317: * Set the context for the time stepper
318: */
319: TSSetApplicationContext(ts, &user);
321: /*
322: * Setup the time stepper, ready for evaluation
323: */
324: TSSetUp(ts);
326: /*
327: * Perform the solve.
328: */
329: TSSolve(ts, x);
330: TSGetSolveTime(ts, &ftime);
331: TSGetStepNumber(ts,&its);
332: PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %D, final time: %4.2e\nResult:\n\n", its, (double)ftime);
333: PrintSolution(x, &user);
335: /*
336: * Free the data structures
337: */
338: VecDestroy(&x);
339: MatDestroy(&J);
340: TSDestroy(&ts);
341: PetscFinalize();
342: return 0;
343: }
345: /*TEST
346: build:
347: requires: !single !complex
349: test:
350: args: -ts_max_steps 8
351: output_file: output/ex42.out
353: TEST*/