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