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