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 { 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(0); 73 } 74 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(0); 138 } 139 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(0); 161 } 162 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(0); 178 } 179 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(0); 192 } 193 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 8 354 output_file: output/ex42.out 355 356 TEST*/ 357