1 static char help[] = "Meinhard't activator-inhibitor model to test TS domain error feature.\n";
2
3 /*
4 The activator-inhibitor on a line is described by the PDE:
5
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
8
9 The PDE part will be solve by finite-difference on the line of cells.
10 */
11
12 #include <petscts.h>
13
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;
25
RHSFunction(TS ts,PetscReal t,Vec X,Vec DXDT,void * ptr)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;
36
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;
47
48 PetscCall(VecGetArrayRead(X, &x));
49 PetscCall(VecGetArray(DXDT, &dxdt));
50
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 }
74
RHSJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void * ptr)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;
86
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;
95
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 }
139
DomainErrorFunction(TS ts,PetscReal t,Vec Y,PetscBool * accept)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;
146
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 }
162
FormInitialState(Vec X,AppCtx * user)163 PetscErrorCode FormInitialState(Vec X, AppCtx *user)
164 {
165 PetscRandom R;
166
167 PetscFunctionBeginUser;
168 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &R));
169 PetscCall(PetscRandomSetFromOptions(R));
170 PetscCall(PetscRandomSetInterval(R, 0., 10.));
171
172 /*
173 * Initialize the state vector
174 */
175 PetscCall(VecSetRandom(X, R));
176 PetscCall(PetscRandomDestroy(&R));
177 PetscFunctionReturn(PETSC_SUCCESS);
178 }
179
PrintSolution(Vec X,AppCtx * user)180 PetscErrorCode PrintSolution(Vec X, AppCtx *user)
181 {
182 const PetscScalar *x;
183 PetscInt i;
184 PetscInt nb_cells = user->nb_cells;
185
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 }
193
main(int argc,char ** argv)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;
203
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");
208
209 /*
210 * Allow user to set the grid dimensions and the equations parameters
211 */
212
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.;
222
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();
234
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));
244
245 /*
246 * Create vector to hold the solution
247 */
248 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2 * user.nb_cells, &x));
249
250 /*
251 * Create time-stepper context
252 */
253 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
254 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
255
256 /*
257 * Tell the time-stepper context where to compute the solution
258 */
259 PetscCall(TSSetSolution(ts, x));
260
261 /*
262 * Allocate the jacobian matrix
263 */
264 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2 * user.nb_cells, 2 * user.nb_cells, 4, 0, &J));
265
266 /*
267 * Provide the call-back for the non-linear function we are evaluating.
268 */
269 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
270
271 /*
272 * Set the Jacobian matrix and the function user to compute Jacobians
273 */
274 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, &user));
275
276 /*
277 * Set the function checking the domain
278 */
279 PetscCall(TSSetFunctionDomainError(ts, DomainErrorFunction));
280
281 /*
282 * Initialize the problem with random values
283 */
284 PetscCall(FormInitialState(x, &user));
285
286 /*
287 * Read the solver type from options
288 */
289 PetscCall(TSSetType(ts, TSPSEUDO));
290
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));
298
299 /*
300 * Set a larger number of potential errors
301 */
302 PetscCall(TSSetMaxStepRejections(ts, 50));
303
304 /*
305 * Also start with a very small dt
306 */
307 PetscCall(TSSetTimeStep(ts, 0.05));
308
309 /*
310 * Set a larger time step increment
311 */
312 PetscCall(TSPseudoSetTimeStepIncrement(ts, 1.5));
313
314 /*
315 * Let the user personalise TS
316 */
317 PetscCall(TSSetFromOptions(ts));
318
319 /*
320 * Set the context for the time stepper
321 */
322 PetscCall(TSSetApplicationContext(ts, &user));
323
324 /*
325 * Setup the time stepper, ready for evaluation
326 */
327 PetscCall(TSSetUp(ts));
328
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));
337
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 }
347
348 /*TEST
349 build:
350 requires: !single !complex
351
352 test:
353 args: -ts_max_steps 20
354 output_file: output/ex42.out
355
356 TEST*/
357