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 PetscErrorCode ierr;
35 PetscScalar *dxdt;
36 const PetscScalar *x;
37
38 PetscFunctionBegin;
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;
48
49 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr);
50 ierr = VecGetArray(DXDT, &dxdt);CHKERRQ(ierr);
51
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 ierr = VecRestoreArray(DXDT, &dxdt);CHKERRQ(ierr);
72 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr);
73 PetscFunctionReturn(0);
74 }
75
RHSJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void * ptr)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;
88
89 PetscFunctionBegin;
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;
97
98 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr);
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 ierr = MatSetValues(B, 1, &rowa, idx, ca, va, INSERT_VALUES);CHKERRQ(ierr);
130 ierr = MatSetValues(B, 1, &rowh, idx, ch, vh, INSERT_VALUES);CHKERRQ(ierr);
131 }
132 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr);
133 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
134 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
135 if (J != B) {
136 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
137 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
138 }
139 PetscFunctionReturn(0);
140 }
141
DomainErrorFunction(TS ts,PetscReal t,Vec Y,PetscBool * accept)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;
149
150 ierr = TSGetApplicationContext(ts, &user);CHKERRQ(ierr);
151 nb_cells = user->nb_cells;
152 ierr = VecGetArrayRead(Y, &x);CHKERRQ(ierr);
153 for (i = 0 ; i < 2*nb_cells ; ++i) {
154 if (PetscRealPart(x[i]) < 0) {
155 ierr = TSGetTimeStep(ts, &dt);CHKERRQ(ierr);
156 ierr = PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t);CHKERRQ(ierr);
157 *accept = PETSC_FALSE;
158 break;
159 }
160 }
161 ierr = VecRestoreArrayRead(Y, &x);CHKERRQ(ierr);
162 PetscFunctionReturn(0);
163 }
164
FormInitialState(Vec X,AppCtx * user)165 PetscErrorCode FormInitialState(Vec X, AppCtx* user)
166 {
167 PetscErrorCode ierr;
168 PetscRandom R;
169
170 PetscFunctionBegin;
171 ierr = PetscRandomCreate(PETSC_COMM_WORLD, &R);CHKERRQ(ierr);
172 ierr = PetscRandomSetFromOptions(R);CHKERRQ(ierr);
173 ierr = PetscRandomSetInterval(R, 0., 10.);CHKERRQ(ierr);
174
175 /*
176 * Initialize the state vector
177 */
178 ierr = VecSetRandom(X, R);CHKERRQ(ierr);
179 ierr = PetscRandomDestroy(&R);CHKERRQ(ierr);
180 PetscFunctionReturn(0);
181 }
182
PrintSolution(Vec X,AppCtx * user)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;
189
190 PetscFunctionBegin;
191 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr);
192 ierr = PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");CHKERRQ(ierr);
193 for (i = 0 ; i < nb_cells ; i++) {
194 ierr = PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2*i], (double)x[2*i+1]);CHKERRQ(ierr);
195 }
196 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr);
197 PetscFunctionReturn(0);
198 }
199
main(int argc,char ** argv)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 */
206 PetscErrorCode ierr;
207 PetscReal ftime;
208 PetscInt its;
209 PetscMPIInt size;
210
211 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
212 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
213 if (size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only");
214
215 /*
216 * Allow user to set the grid dimensions and the equations parameters
217 */
218
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.;
228
229 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem settings", "PROBLEM");CHKERRQ(ierr);
230 ierr = PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c",user.nb_cells, &user.nb_cells,NULL);CHKERRQ(ierr);
231 ierr = PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c",user.alpha, &user.alpha,NULL);CHKERRQ(ierr);
232 ierr = PetscOptionsReal("-beta", "Inhibition factor", "ex42.c",user.beta, &user.beta,NULL);CHKERRQ(ierr);
233 ierr = PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c",user.rho_a, &user.rho_a,NULL);CHKERRQ(ierr);
234 ierr = PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c",user.mu_a, &user.mu_a,NULL);CHKERRQ(ierr);
235 ierr = PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c",user.D_a, &user.D_a,NULL);CHKERRQ(ierr);
236 ierr = PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c",user.rho_h, &user.rho_h,NULL);CHKERRQ(ierr);
237 ierr = PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c",user.mu_h, &user.mu_h,NULL);CHKERRQ(ierr);
238 ierr = PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c",user.D_h, &user.D_h,NULL);CHKERRQ(ierr);
239 ierr = PetscOptionsEnd();CHKERRQ(ierr);
240
241 ierr = PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %D\n", user.nb_cells);CHKERRQ(ierr);
242 ierr = PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", (double)user.alpha);CHKERRQ(ierr);
243 ierr = PetscPrintf(PETSC_COMM_WORLD, "beta: %5.5g\n", (double)user.beta);CHKERRQ(ierr);
244 ierr = PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", (double)user.rho_a);CHKERRQ(ierr);
245 ierr = PetscPrintf(PETSC_COMM_WORLD, "mu_a: %5.5g\n", (double)user.mu_a);CHKERRQ(ierr);
246 ierr = PetscPrintf(PETSC_COMM_WORLD, "D_a: %5.5g\n", (double)user.D_a);CHKERRQ(ierr);
247 ierr = PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", (double)user.rho_h);CHKERRQ(ierr);
248 ierr = PetscPrintf(PETSC_COMM_WORLD, "mu_h: %5.5g\n", (double)user.mu_h);CHKERRQ(ierr);
249 ierr = PetscPrintf(PETSC_COMM_WORLD, "D_h: %5.5g\n", (double)user.D_h);CHKERRQ(ierr);
250
251 /*
252 * Create vector to hold the solution
253 */
254 ierr = VecCreateSeq(PETSC_COMM_WORLD, 2*user.nb_cells, &x);CHKERRQ(ierr);
255
256 /*
257 * Create time-stepper context
258 */
259 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
260 ierr = TSSetProblemType(ts, TS_NONLINEAR);CHKERRQ(ierr);
261
262 /*
263 * Tell the time-stepper context where to compute the solution
264 */
265 ierr = TSSetSolution(ts, x);CHKERRQ(ierr);
266
267 /*
268 * Allocate the jacobian matrix
269 */
270 ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD, 2*user.nb_cells, 2*user.nb_cells, 4, 0, &J);CHKERRQ(ierr);
271
272 /*
273 * Provide the call-back for the non-linear function we are evaluating.
274 */
275 ierr = TSSetRHSFunction(ts, NULL, RHSFunction, &user);CHKERRQ(ierr);
276
277 /*
278 * Set the Jacobian matrix and the function user to compute Jacobians
279 */
280 ierr = TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);CHKERRQ(ierr);
281
282 /*
283 * Set the function checking the domain
284 */
285 ierr = TSSetFunctionDomainError(ts, &DomainErrorFunction);CHKERRQ(ierr);
286
287 /*
288 * Initialize the problem with random values
289 */
290 ierr = FormInitialState(x, &user);CHKERRQ(ierr);
291
292 /*
293 * Read the solver type from options
294 */
295 ierr = TSSetType(ts, TSPSEUDO);CHKERRQ(ierr);
296
297 /*
298 * Set a large number of timesteps and final duration time to insure
299 * convergenge to steady state
300 */
301 ierr = TSSetMaxSteps(ts, 2147483647);CHKERRQ(ierr);
302 ierr = TSSetMaxTime(ts, 1.e12);CHKERRQ(ierr);
303 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
304
305 /*
306 * Set a larger number of potential errors
307 */
308 ierr = TSSetMaxStepRejections(ts, 50);CHKERRQ(ierr);
309
310 /*
311 * Also start with a very small dt
312 */
313 ierr = TSSetTimeStep(ts, 0.05);CHKERRQ(ierr);
314
315 /*
316 * Set a larger time step increment
317 */
318 ierr = TSPseudoSetTimeStepIncrement(ts, 1.5);CHKERRQ(ierr);
319
320 /*
321 * Let the user personalise TS
322 */
323 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
324
325 /*
326 * Set the context for the time stepper
327 */
328 ierr = TSSetApplicationContext(ts, &user);CHKERRQ(ierr);
329
330 /*
331 * Setup the time stepper, ready for evaluation
332 */
333 ierr = TSSetUp(ts);CHKERRQ(ierr);
334
335 /*
336 * Perform the solve.
337 */
338 ierr = TSSolve(ts, x);CHKERRQ(ierr);
339 ierr = TSGetSolveTime(ts, &ftime);CHKERRQ(ierr);
340 ierr = TSGetStepNumber(ts,&its);CHKERRQ(ierr);
341 ierr = PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %D, final time: %4.2e\nResult:\n\n", its, (double)ftime);CHKERRQ(ierr);
342 ierr = PrintSolution(x, &user);CHKERRQ(ierr);
343
344 /*
345 * Free the data structures
346 */
347 ierr = VecDestroy(&x);CHKERRQ(ierr);
348 ierr = MatDestroy(&J);CHKERRQ(ierr);
349 ierr = TSDestroy(&ts);CHKERRQ(ierr);
350 ierr = PetscFinalize();
351 return ierr;
352 }
353
354 /*TEST
355 build:
356 requires: !single !complex
357
358 test:
359 args: -ts_max_steps 8
360 output_file: output/ex42.out
361
362 TEST*/
363