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: PetscFunctionBeginUser;
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: PetscCall(VecGetArrayRead(X, &x));
49: PetscCall(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: PetscCall(VecRestoreArray(DXDT, &dxdt));
71: PetscCall(VecRestoreArrayRead(X, &x));
72: PetscFunctionReturn(PETSC_SUCCESS);
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;
87: PetscFunctionBeginUser;
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: PetscCall(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: PetscCall(MatSetValues(B, 1, &rowa, idx, ca, va, INSERT_VALUES));
128: PetscCall(MatSetValues(B, 1, &rowh, idx, ch, vh, INSERT_VALUES));
129: }
130: PetscCall(VecRestoreArrayRead(X, &x));
131: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
132: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
133: if (J != B) {
134: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
135: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
136: }
137: PetscFunctionReturn(PETSC_SUCCESS);
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;
147: PetscFunctionBeginUser;
148: PetscCall(TSGetApplicationContext(ts, &user));
149: nb_cells = user->nb_cells;
150: PetscCall(VecGetArrayRead(Y, &x));
151: for (i = 0; i < 2 * nb_cells; ++i) {
152: if (PetscRealPart(x[i]) < 0) {
153: PetscCall(TSGetTimeStep(ts, &dt));
154: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t));
155: *accept = PETSC_FALSE;
156: break;
157: }
158: }
159: PetscCall(VecRestoreArrayRead(Y, &x));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: PetscErrorCode FormInitialState(Vec X, AppCtx *user)
164: {
165: PetscRandom R;
167: PetscFunctionBeginUser;
168: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &R));
169: PetscCall(PetscRandomSetFromOptions(R));
170: PetscCall(PetscRandomSetInterval(R, 0., 10.));
172: /*
173: * Initialize the state vector
174: */
175: PetscCall(VecSetRandom(X, R));
176: PetscCall(PetscRandomDestroy(&R));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
181: {
182: const PetscScalar *x;
183: PetscInt i;
184: PetscInt nb_cells = user->nb_cells;
186: PetscFunctionBeginUser;
187: PetscCall(VecGetArrayRead(X, &x));
188: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n"));
189: for (i = 0; i < nb_cells; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2 * i], (double)x[2 * i + 1]));
190: PetscCall(VecRestoreArrayRead(X, &x));
191: PetscFunctionReturn(PETSC_SUCCESS);
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;
204: PetscFunctionBeginUser;
205: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
206: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
207: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
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: PetscCall(PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c", user.nb_cells, &user.nb_cells, NULL));
225: PetscCall(PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c", user.alpha, &user.alpha, NULL));
226: PetscCall(PetscOptionsReal("-beta", "Inhibition factor", "ex42.c", user.beta, &user.beta, NULL));
227: PetscCall(PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c", user.rho_a, &user.rho_a, NULL));
228: PetscCall(PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c", user.mu_a, &user.mu_a, NULL));
229: PetscCall(PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c", user.D_a, &user.D_a, NULL));
230: PetscCall(PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c", user.rho_h, &user.rho_h, NULL));
231: PetscCall(PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c", user.mu_h, &user.mu_h, NULL));
232: PetscCall(PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c", user.D_h, &user.D_h, NULL));
233: PetscOptionsEnd();
235: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %" PetscInt_FMT "\n", user.nb_cells));
236: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", (double)user.alpha));
237: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "beta: %5.5g\n", (double)user.beta));
238: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", (double)user.rho_a));
239: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu_a: %5.5g\n", (double)user.mu_a));
240: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "D_a: %5.5g\n", (double)user.D_a));
241: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", (double)user.rho_h));
242: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu_h: %5.5g\n", (double)user.mu_h));
243: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "D_h: %5.5g\n", (double)user.D_h));
245: /*
246: * Create vector to hold the solution
247: */
248: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2 * user.nb_cells, &x));
250: /*
251: * Create time-stepper context
252: */
253: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
254: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
256: /*
257: * Tell the time-stepper context where to compute the solution
258: */
259: PetscCall(TSSetSolution(ts, x));
261: /*
262: * Allocate the jacobian matrix
263: */
264: PetscCall(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: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
271: /*
272: * Set the Jacobian matrix and the function user to compute Jacobians
273: */
274: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, &user));
276: /*
277: * Set the function checking the domain
278: */
279: PetscCall(TSSetFunctionDomainError(ts, &DomainErrorFunction));
281: /*
282: * Initialize the problem with random values
283: */
284: PetscCall(FormInitialState(x, &user));
286: /*
287: * Read the solver type from options
288: */
289: PetscCall(TSSetType(ts, TSPSEUDO));
291: /*
292: * Set a large number of timesteps and final duration time to insure
293: * convergenge to steady state
294: */
295: PetscCall(TSSetMaxSteps(ts, 2147483647));
296: PetscCall(TSSetMaxTime(ts, 1.e12));
297: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
299: /*
300: * Set a larger number of potential errors
301: */
302: PetscCall(TSSetMaxStepRejections(ts, 50));
304: /*
305: * Also start with a very small dt
306: */
307: PetscCall(TSSetTimeStep(ts, 0.05));
309: /*
310: * Set a larger time step increment
311: */
312: PetscCall(TSPseudoSetTimeStepIncrement(ts, 1.5));
314: /*
315: * Let the user personalise TS
316: */
317: PetscCall(TSSetFromOptions(ts));
319: /*
320: * Set the context for the time stepper
321: */
322: PetscCall(TSSetApplicationContext(ts, &user));
324: /*
325: * Setup the time stepper, ready for evaluation
326: */
327: PetscCall(TSSetUp(ts));
329: /*
330: * Perform the solve.
331: */
332: PetscCall(TSSolve(ts, x));
333: PetscCall(TSGetSolveTime(ts, &ftime));
334: PetscCall(TSGetStepNumber(ts, &its));
335: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %" PetscInt_FMT ", final time: %4.2e\nResult:\n\n", its, (double)ftime));
336: PetscCall(PrintSolution(x, &user));
338: /*
339: * Free the data structures
340: */
341: PetscCall(VecDestroy(&x));
342: PetscCall(MatDestroy(&J));
343: PetscCall(TSDestroy(&ts));
344: PetscCall(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*/