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: 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*/