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;
151: TSGetApplicationContext(ts, &user);
152: nb_cells = user->nb_cells;
153: VecGetArrayRead(Y, &x);
154: for (i = 0 ; i < 2*nb_cells ; ++i) {
155: if (PetscRealPart(x[i]) < 0) {
156: TSGetTimeStep(ts, &dt);
157: PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t);
158: *accept = PETSC_FALSE;
159: break;
160: }
161: }
162: VecRestoreArrayRead(Y, &x);
163: return(0);
164: }
166: PetscErrorCode FormInitialState(Vec X, AppCtx* user)
167: {
169: PetscRandom R;
172: PetscRandomCreate(PETSC_COMM_WORLD, &R);
173: PetscRandomSetFromOptions(R);
174: PetscRandomSetInterval(R, 0., 10.);
176: /*
177: * Initialize the state vector
178: */
179: VecSetRandom(X, R);
180: PetscRandomDestroy(&R);
181: return(0);
182: }
184: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
185: {
186: PetscErrorCode ierr;
187: const PetscScalar *x;
188: PetscInt i;
189: PetscInt nb_cells = user->nb_cells;
192: VecGetArrayRead(X, &x);
193: PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");
194: for (i = 0 ; i < nb_cells ; i++) {
195: PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2*i], (double)x[2*i+1]);
196: }
197: VecRestoreArrayRead(X, &x);
198: return(0);
199: }
201: int main(int argc, char **argv)
202: {
203: TS ts; /* time-stepping context */
204: Vec x; /* State vector */
205: Mat J; /* Jacobian matrix */
206: AppCtx user; /* user-defined context */
208: PetscReal ftime;
209: PetscInt its;
210: PetscMPIInt size;
212: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
213: MPI_Comm_size(PETSC_COMM_WORLD, &size);
214: if (size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only");
216: /*
217: * Allow user to set the grid dimensions and the equations parameters
218: */
220: user.nb_cells = 50;
221: user.alpha = 10.;
222: user.beta = 1.;
223: user.rho_a = 1.;
224: user.rho_h = 2.;
225: user.mu_a = 2.;
226: user.mu_h = 3.;
227: user.D_a = 0.;
228: user.D_h = 30.;
230: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem settings", "PROBLEM");
231: PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c",user.nb_cells, &user.nb_cells,NULL);
232: PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c",user.alpha, &user.alpha,NULL);
233: PetscOptionsReal("-beta", "Inhibition factor", "ex42.c",user.beta, &user.beta,NULL);
234: PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c",user.rho_a, &user.rho_a,NULL);
235: PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c",user.mu_a, &user.mu_a,NULL);
236: PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c",user.D_a, &user.D_a,NULL);
237: PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c",user.rho_h, &user.rho_h,NULL);
238: PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c",user.mu_h, &user.mu_h,NULL);
239: PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c",user.D_h, &user.D_h,NULL);
240: PetscOptionsEnd();
242: PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %D\n", user.nb_cells);
243: PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", (double)user.alpha);
244: PetscPrintf(PETSC_COMM_WORLD, "beta: %5.5g\n", (double)user.beta);
245: PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", (double)user.rho_a);
246: PetscPrintf(PETSC_COMM_WORLD, "mu_a: %5.5g\n", (double)user.mu_a);
247: PetscPrintf(PETSC_COMM_WORLD, "D_a: %5.5g\n", (double)user.D_a);
248: PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", (double)user.rho_h);
249: PetscPrintf(PETSC_COMM_WORLD, "mu_h: %5.5g\n", (double)user.mu_h);
250: PetscPrintf(PETSC_COMM_WORLD, "D_h: %5.5g\n", (double)user.D_h);
252: /*
253: * Create vector to hold the solution
254: */
255: VecCreateSeq(PETSC_COMM_WORLD, 2*user.nb_cells, &x);
257: /*
258: * Create time-stepper context
259: */
260: TSCreate(PETSC_COMM_WORLD, &ts);
261: TSSetProblemType(ts, TS_NONLINEAR);
263: /*
264: * Tell the time-stepper context where to compute the solution
265: */
266: TSSetSolution(ts, x);
268: /*
269: * Allocate the jacobian matrix
270: */
271: MatCreateSeqAIJ(PETSC_COMM_WORLD, 2*user.nb_cells, 2*user.nb_cells, 4, 0, &J);
273: /*
274: * Provide the call-back for the non-linear function we are evaluating.
275: */
276: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
278: /*
279: * Set the Jacobian matrix and the function user to compute Jacobians
280: */
281: TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);
283: /*
284: * Set the function checking the domain
285: */
286: TSSetFunctionDomainError(ts, &DomainErrorFunction);
288: /*
289: * Initialize the problem with random values
290: */
291: FormInitialState(x, &user);
293: /*
294: * Read the solver type from options
295: */
296: TSSetType(ts, TSPSEUDO);
298: /*
299: * Set a large number of timesteps and final duration time to insure
300: * convergenge to steady state
301: */
302: TSSetMaxSteps(ts, 2147483647);
303: TSSetMaxTime(ts, 1.e12);
304: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
306: /*
307: * Set a larger number of potential errors
308: */
309: TSSetMaxStepRejections(ts, 50);
311: /*
312: * Also start with a very small dt
313: */
314: TSSetTimeStep(ts, 0.05);
316: /*
317: * Set a larger time step increment
318: */
319: TSPseudoSetTimeStepIncrement(ts, 1.5);
321: /*
322: * Let the user personalise TS
323: */
324: TSSetFromOptions(ts);
326: /*
327: * Set the context for the time stepper
328: */
329: TSSetApplicationContext(ts, &user);
331: /*
332: * Setup the time stepper, ready for evaluation
333: */
334: TSSetUp(ts);
336: /*
337: * Perform the solve.
338: */
339: TSSolve(ts, x);
340: TSGetSolveTime(ts, &ftime);
341: TSGetStepNumber(ts,&its);
342: PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %D, final time: %4.2e\nResult:\n\n", its, (double)ftime);
343: PrintSolution(x, &user);
345: /*
346: * Free the data structures
347: */
348: VecDestroy(&x);
349: MatDestroy(&J);
350: TSDestroy(&ts);
351: PetscFinalize();
352: return ierr;
353: }
355: /*TEST
356: build:
357: requires: !single !complex
359: test:
360: args: -ts_max_steps 8
361: output_file: output/ex42.out
363: TEST*/