Actual source code: ex42.c
petsc-3.7.5 2017-01-01
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;
28: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec DXDT, void* ptr)
29: {
30: AppCtx* user = (AppCtx*)ptr;
31: PetscInt nb_cells, i;
32: PetscReal alpha, beta;
33: PetscReal rho_a, mu_a, D_a;
34: PetscReal rho_h, mu_h, D_h;
35: PetscReal a, h, da, dh, d2a, d2h;
36: PetscErrorCode ierr;
37: PetscScalar *dxdt;
38: const PetscScalar *x;
41: nb_cells = user->nb_cells;
42: alpha = user->alpha;
43: beta = user->beta;
44: rho_a = user->rho_a;
45: mu_a = user->mu_a;
46: D_a = user->D_a;
47: rho_h = user->rho_h;
48: mu_h = user->mu_h;
49: D_h = user->D_h;
51: VecGetArrayRead(X, &x);
52: VecGetArray(DXDT, &dxdt);
54: for(i = 0 ; i < nb_cells ; i++) {
55: a = x[2*i];
56: h = x[2*i+1];
57: // Reaction:
58: da = alpha * a*a / (1. + beta * h) + rho_a - mu_a * a;
59: dh = alpha * a*a + rho_h - mu_h*h;
60: // Diffusion:
61: d2a = d2h = 0.;
62: if(i > 0) {
63: d2a += (x[2*(i-1)] - a);
64: d2h += (x[2*(i-1)+1] - h);
65: }
66: if(i < nb_cells-1) {
67: d2a += (x[2*(i+1)] - a);
68: d2h += (x[2*(i+1)+1] - h);
69: }
70: dxdt[2*i] = da + D_a*d2a;
71: dxdt[2*i+1] = dh + D_h*d2h;
72: }
73: VecRestoreArray(DXDT, &dxdt);
74: VecRestoreArrayRead(X, &x);
75: return(0);
76: }
80: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
81: {
82: AppCtx *user = (AppCtx*)ptr;
83: PetscInt nb_cells, i, idx;
84: PetscReal alpha, beta;
85: PetscReal mu_a, D_a;
86: PetscReal mu_h, D_h;
87: PetscReal a, h;
88: const PetscScalar *x;
89: PetscScalar va[4], vh[4];
90: PetscInt ca[4], ch[4], rowa, rowh;
91: PetscErrorCode ierr;
94: nb_cells = user->nb_cells;
95: alpha = user->alpha;
96: beta = user->beta;
97: mu_a = user->mu_a;
98: D_a = user->D_a;
99: mu_h = user->mu_h;
100: D_h = user->D_h;
102: VecGetArrayRead(X, &x);
103: for(i = 0; i < nb_cells ; ++i) {
104: rowa = 2*i;
105: rowh = 2*i+1;
106: a = x[2*i];
107: h = x[2*i+1];
108: ca[0] = ch[1] = 2*i;
109: va[0] = 2*alpha*a / (1.+beta*h) - mu_a;
110: vh[1] = 2*alpha*a;
111: ca[1] = ch[0] = 2*i+1;
112: va[1] = -beta*alpha*a*a / ((1.+beta*h)*(1.+beta*h));
113: vh[0] = -mu_h;
114: idx = 2;
115: if(i > 0) {
116: ca[idx] = 2*(i-1);
117: ch[idx] = 2*(i-1)+1;
118: va[idx] = D_a;
119: vh[idx] = D_h;
120: va[0] -= D_a;
121: vh[0] -= D_h;
122: idx++;
123: }
124: if(i < nb_cells-1) {
125: ca[idx] = 2*(i+1);
126: ch[idx] = 2*(i+1)+1;
127: va[idx] = D_a;
128: vh[idx] = D_h;
129: va[0] -= D_a;
130: vh[0] -= D_h;
131: idx++;
132: }
133: MatSetValues(B, 1, &rowa, idx, ca, va, INSERT_VALUES);
134: MatSetValues(B, 1, &rowh, idx, ch, vh, INSERT_VALUES);
135: }
136: VecRestoreArrayRead(X, &x);
137: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
139: if (J != B) {
140: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
141: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
142: }
143: return(0);
144: }
148: PetscErrorCode DomainErrorFunction(TS ts, PetscReal t, Vec Y, PetscBool *accept)
149: {
150: AppCtx *user;
151: PetscReal dt;
152: PetscErrorCode ierr;
153: const PetscScalar *x;
154: PetscInt nb_cells, i;
156: TSGetApplicationContext(ts, &user);
157: nb_cells = user->nb_cells;
158: VecGetArrayRead(Y, &x);
159: for(i = 0 ; i < 2*nb_cells ; ++i) {
160: if(PetscRealPart(x[i]) < 0) {
161: TSGetTimeStep(ts, &dt);
162: PetscPrintf(PETSC_COMM_WORLD, " ** Domain Error at time %g\n", (double)t);
163: *accept = PETSC_FALSE;
164: break;
165: }
166: }
167: VecRestoreArrayRead(Y, &x);
168: return(0);
169: }
173: PetscErrorCode FormInitialState(Vec X, AppCtx* user)
174: {
176: PetscRandom R;
179: PetscRandomCreate(PETSC_COMM_WORLD, &R);
180: PetscRandomSetFromOptions(R);
181: PetscRandomSetInterval(R, 0., 10.);
183: /*
184: * Initialize the state vector
185: */
186: VecSetRandom(X, R);
187: PetscRandomDestroy(&R);
188: return(0);
189: }
193: PetscErrorCode PrintSolution(Vec X, AppCtx *user)
194: {
195: PetscErrorCode ierr;
196: const PetscScalar *x;
197: PetscInt i;
198: PetscInt nb_cells = user->nb_cells;
201: VecGetArrayRead(X, &x);
202: PetscPrintf(PETSC_COMM_WORLD, "Activator,Inhibitor\n");
203: for(i = 0 ; i < nb_cells ; i++) {
204: PetscPrintf(PETSC_COMM_WORLD, "%5.6e,%5.6e\n", (double)x[2*i], (double)x[2*i+1]);
205: }
206: VecRestoreArrayRead(X, &x);
207: return(0);
208: }
212: int main(int argc, char **argv)
213: {
214: TS ts; /* time-stepping context */
215: Vec x; /* State vector */
216: Mat J; /* Jacobian matrix */
217: AppCtx user; /* user-defined context */
219: PetscReal ftime;
220: PetscInt its;
221: PetscMPIInt size;
223: PetscInitialize(&argc, &argv, NULL, help);
224: MPI_Comm_size(PETSC_COMM_WORLD, &size);
225: if(size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only");
227: /*
228: * Allow user to set the grid dimensions and the equations parameters
229: */
231: user.nb_cells = 50;
232: user.alpha = 10.;
233: user.beta = 1.;
234: user.rho_a = 1.;
235: user.rho_h = 2.;
236: user.mu_a = 2.;
237: user.mu_h = 3.;
238: user.D_a = 0.;
239: user.D_h = 30.;
241: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem settings", "PROBLEM");
242: PetscOptionsInt("-nb_cells", "Number of cells", "ex42.c",user.nb_cells, &user.nb_cells,NULL);
243: PetscOptionsReal("-alpha", "Autocatalysis factor", "ex42.c",user.alpha, &user.alpha,NULL);
244: PetscOptionsReal("-beta", "Inhibition factor", "ex42.c",user.beta, &user.beta,NULL);
245: PetscOptionsReal("-rho_a", "Default production of the activator", "ex42.c",user.rho_a, &user.rho_a,NULL);
246: PetscOptionsReal("-mu_a", "Degradation rate of the activator", "ex42.c",user.mu_a, &user.mu_a,NULL);
247: PetscOptionsReal("-D_a", "Diffusion rate of the activator", "ex42.c",user.D_a, &user.D_a,NULL);
248: PetscOptionsReal("-rho_h", "Default production of the inhibitor", "ex42.c",user.rho_h, &user.rho_h,NULL);
249: PetscOptionsReal("-mu_h", "Degradation rate of the inhibitor", "ex42.c",user.mu_h, &user.mu_h,NULL);
250: PetscOptionsReal("-D_h", "Diffusion rate of the inhibitor", "ex42.c",user.D_h, &user.D_h,NULL);
251: PetscOptionsEnd();
253: PetscPrintf(PETSC_COMM_WORLD, "nb_cells: %D\n", user.nb_cells);
254: PetscPrintf(PETSC_COMM_WORLD, "alpha: %5.5g\n", user.alpha);
255: PetscPrintf(PETSC_COMM_WORLD, "beta: %5.5g\n", user.beta);
256: PetscPrintf(PETSC_COMM_WORLD, "rho_a: %5.5g\n", user.rho_a);
257: PetscPrintf(PETSC_COMM_WORLD, "mu_a: %5.5g\n", user.mu_a);
258: PetscPrintf(PETSC_COMM_WORLD, "D_a: %5.5g\n", user.D_a);
259: PetscPrintf(PETSC_COMM_WORLD, "rho_h: %5.5g\n", user.rho_h);
260: PetscPrintf(PETSC_COMM_WORLD, "mu_h: %5.5g\n", user.mu_h);
261: PetscPrintf(PETSC_COMM_WORLD, "D_h: %5.5g\n", user.D_h);
263: /*
264: * Create vector to hold the solution
265: */
266: VecCreateSeq(PETSC_COMM_WORLD, 2*user.nb_cells, &x);
268: /*
269: * Create time-stepper context
270: */
271: TSCreate(PETSC_COMM_WORLD, &ts);
272: TSSetProblemType(ts, TS_NONLINEAR);
274: /*
275: * Tell the time-stepper context where to compute the solution
276: */
277: TSSetSolution(ts, x);
279: /*
280: * Allocate the jacobian matrix
281: */
282: MatCreateSeqAIJ(PETSC_COMM_WORLD, 2*user.nb_cells, 2*user.nb_cells, 4, 0, &J);
284: /*
285: * Provide the call-back for the non-linear function we are evaluating.
286: */
287: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
289: /*
290: * Set the Jacobian matrix and the function user to compute Jacobians
291: */
292: TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);
294: /*
295: * Set the function checking the domain
296: */
297: TSSetFunctionDomainError(ts, &DomainErrorFunction);
299: /*
300: * Initialize the problem with random values
301: */
302: FormInitialState(x, &user);
304: /*
305: * Read the solver type from options
306: */
307: TSSetType(ts, TSPSEUDO);
309: /*
310: * Set a large number of timesteps and final duration time to insure
311: * convergenge to steady state
312: */
313: TSSetDuration(ts, 5000, 1e12);
314: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
316: /*
317: * Set a larger number of potential errors
318: */
319: TSSetMaxStepRejections(ts, 50);
321: /*
322: * Also start with a very small dt
323: */
324: TSSetTimeStep(ts, 0.05);
326: /*
327: * Set a larger time step increment
328: */
329: TSPseudoSetTimeStepIncrement(ts, 1.5);
331: /*
332: * Let the user personalise TS
333: */
334: TSSetFromOptions(ts);
336: /*
337: * Set the context for the time stepper
338: */
339: TSSetApplicationContext(ts, &user);
341: /*
342: * Setup the time stepper, ready for evaluation
343: */
344: TSSetUp(ts);
346: /*
347: * Perform the solve.
348: */
349: TSSolve(ts, x);
350: TSGetSolveTime(ts, &ftime);
351: TSGetTimeStepNumber(ts,&its);
352: PetscPrintf(PETSC_COMM_WORLD, "Number of time steps = %D, final time: %4.2e\nResult:\n\n", its, (double)ftime);
353: PrintSolution(x, &user);
355: /*
356: * Free the data structures
357: */
358: VecDestroy(&x);
359: MatDestroy(&J);
360: TSDestroy(&ts);
361: PetscFinalize();
362: return 0;
363: }