1 static const char help[] = "Time-dependent advection-reaction PDE in 1d, demonstrates IMEX methods.\n"; 2 /* 3 u_t + a1*u_x = -k1*u + k2*v + s1 4 v_t + a2*v_x = k1*u - k2*v + s2 5 0 < x < 1; 6 a1 = 1, k1 = 10^6, s1 = 0, 7 a2 = 0, k2 = 2*k1, s2 = 1 8 9 Initial conditions: 10 u(x,0) = 1 + s2*x 11 v(x,0) = k0/k1*u(x,0) + s1/k1 12 13 Upstream boundary conditions: 14 u(0,t) = 1-sin(12*t)^4 15 16 Useful command-line parameters: 17 -ts_arkimex_fully_implicit - treats advection implicitly, only relevant with -ts_type arkimex (default) 18 -ts_type rosw - use Rosenbrock-W method (linearly implicit IMEX, amortizes assembly and preconditioner setup) 19 */ 20 21 #include <petscdm.h> 22 #include <petscdmda.h> 23 #include <petscts.h> 24 25 typedef PetscScalar Field[2]; 26 27 typedef struct _User *User; 28 struct _User { 29 PetscReal a[2]; /* Advection speeds */ 30 PetscReal k[2]; /* Reaction coefficients */ 31 PetscReal s[2]; /* Source terms */ 32 }; 33 34 static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *); 35 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *); 36 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 37 static PetscErrorCode FormInitialSolution(TS, Vec, void *); 38 39 int main(int argc, char **argv) { 40 TS ts; /* time integrator */ 41 SNES snes; /* nonlinear solver */ 42 SNESLineSearch linesearch; /* line search */ 43 Vec X; /* solution, residual vectors */ 44 Mat J; /* Jacobian matrix */ 45 PetscInt steps, mx; 46 DM da; 47 PetscReal ftime, dt; 48 struct _User user; /* user-defined work context */ 49 TSConvergedReason reason; 50 51 PetscFunctionBeginUser; 52 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 53 54 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 55 Create distributed array (DMDA) to manage parallel grid and vectors 56 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 57 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da)); 58 PetscCall(DMSetFromOptions(da)); 59 PetscCall(DMSetUp(da)); 60 61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 Extract global vectors from DMDA; 63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 64 PetscCall(DMCreateGlobalVector(da, &X)); 65 66 /* Initialize user application context */ 67 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", ""); 68 { 69 user.a[0] = 1; 70 PetscCall(PetscOptionsReal("-a0", "Advection rate 0", "", user.a[0], &user.a[0], NULL)); 71 user.a[1] = 0; 72 PetscCall(PetscOptionsReal("-a1", "Advection rate 1", "", user.a[1], &user.a[1], NULL)); 73 user.k[0] = 1e6; 74 PetscCall(PetscOptionsReal("-k0", "Reaction rate 0", "", user.k[0], &user.k[0], NULL)); 75 user.k[1] = 2 * user.k[0]; 76 PetscCall(PetscOptionsReal("-k1", "Reaction rate 1", "", user.k[1], &user.k[1], NULL)); 77 user.s[0] = 0; 78 PetscCall(PetscOptionsReal("-s0", "Source 0", "", user.s[0], &user.s[0], NULL)); 79 user.s[1] = 1; 80 PetscCall(PetscOptionsReal("-s1", "Source 1", "", user.s[1], &user.s[1], NULL)); 81 } 82 PetscOptionsEnd(); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Create timestepping solver context 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 88 PetscCall(TSSetDM(ts, da)); 89 PetscCall(TSSetType(ts, TSARKIMEX)); 90 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user)); 91 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user)); 92 PetscCall(DMSetMatType(da, MATAIJ)); 93 PetscCall(DMCreateMatrix(da, &J)); 94 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user)); 95 96 /* A line search in the nonlinear solve can fail due to ill-conditioning unless an absolute tolerance is set. Since 97 * this problem is linear, we deactivate the line search. For a linear problem, it is usually recommended to also use 98 * SNESSetType(snes,SNESKSPONLY). */ 99 PetscCall(TSGetSNES(ts, &snes)); 100 PetscCall(SNESGetLineSearch(snes, &linesearch)); 101 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC)); 102 103 ftime = .1; 104 PetscCall(TSSetMaxTime(ts, ftime)); 105 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 106 107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 108 Set initial conditions 109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 110 PetscCall(FormInitialSolution(ts, X, &user)); 111 PetscCall(TSSetSolution(ts, X)); 112 PetscCall(VecGetSize(X, &mx)); 113 dt = .1 * PetscMax(user.a[0], user.a[1]) / mx; /* Advective CFL, I don't know why it needs so much safety factor. */ 114 PetscCall(TSSetTimeStep(ts, dt)); 115 116 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 117 Set runtime options 118 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 119 PetscCall(TSSetFromOptions(ts)); 120 121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122 Solve nonlinear system 123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 124 PetscCall(TSSolve(ts, X)); 125 PetscCall(TSGetSolveTime(ts, &ftime)); 126 PetscCall(TSGetStepNumber(ts, &steps)); 127 PetscCall(TSGetConvergedReason(ts, &reason)); 128 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps)); 129 130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 Free work space. 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 PetscCall(MatDestroy(&J)); 134 PetscCall(VecDestroy(&X)); 135 PetscCall(TSDestroy(&ts)); 136 PetscCall(DMDestroy(&da)); 137 PetscCall(PetscFinalize()); 138 return 0; 139 } 140 141 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr) { 142 User user = (User)ptr; 143 DM da; 144 DMDALocalInfo info; 145 PetscInt i; 146 Field *f; 147 const Field *x, *xdot; 148 149 PetscFunctionBeginUser; 150 PetscCall(TSGetDM(ts, &da)); 151 PetscCall(DMDAGetLocalInfo(da, &info)); 152 153 /* Get pointers to vector data */ 154 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x)); 155 PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot)); 156 PetscCall(DMDAVecGetArray(da, F, &f)); 157 158 /* Compute function over the locally owned part of the grid */ 159 for (i = info.xs; i < info.xs + info.xm; i++) { 160 f[i][0] = xdot[i][0] + user->k[0] * x[i][0] - user->k[1] * x[i][1] - user->s[0]; 161 f[i][1] = xdot[i][1] - user->k[0] * x[i][0] + user->k[1] * x[i][1] - user->s[1]; 162 } 163 164 /* Restore vectors */ 165 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x)); 166 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot)); 167 PetscCall(DMDAVecRestoreArray(da, F, &f)); 168 PetscFunctionReturn(0); 169 } 170 171 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 172 User user = (User)ptr; 173 DM da; 174 Vec Xloc; 175 DMDALocalInfo info; 176 PetscInt i, j; 177 PetscReal hx; 178 Field *f; 179 const Field *x; 180 181 PetscFunctionBeginUser; 182 PetscCall(TSGetDM(ts, &da)); 183 PetscCall(DMDAGetLocalInfo(da, &info)); 184 hx = 1.0 / (PetscReal)info.mx; 185 186 /* 187 Scatter ghost points to local vector,using the 2-step process 188 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 189 By placing code between these two statements, computations can be 190 done while messages are in transition. 191 */ 192 PetscCall(DMGetLocalVector(da, &Xloc)); 193 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 194 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 195 196 /* Get pointers to vector data */ 197 PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x)); 198 PetscCall(DMDAVecGetArray(da, F, &f)); 199 200 /* Compute function over the locally owned part of the grid */ 201 for (i = info.xs; i < info.xs + info.xm; i++) { 202 const PetscReal *a = user->a; 203 PetscReal u0t[2]; 204 u0t[0] = 1.0 - PetscPowRealInt(PetscSinReal(12 * t), 4); 205 u0t[1] = 0.0; 206 for (j = 0; j < 2; j++) { 207 if (i == 0) f[i][j] = a[j] / hx * (1. / 3 * u0t[j] + 0.5 * x[i][j] - x[i + 1][j] + 1. / 6 * x[i + 2][j]); 208 else if (i == 1) f[i][j] = a[j] / hx * (-1. / 12 * u0t[j] + 2. / 3 * x[i - 1][j] - 2. / 3 * x[i + 1][j] + 1. / 12 * x[i + 2][j]); 209 else if (i == info.mx - 2) f[i][j] = a[j] / hx * (-1. / 6 * x[i - 2][j] + x[i - 1][j] - 0.5 * x[i][j] - 1. / 3 * x[i + 1][j]); 210 else if (i == info.mx - 1) f[i][j] = a[j] / hx * (-x[i][j] + x[i - 1][j]); 211 else f[i][j] = a[j] / hx * (-1. / 12 * x[i - 2][j] + 2. / 3 * x[i - 1][j] - 2. / 3 * x[i + 1][j] + 1. / 12 * x[i + 2][j]); 212 } 213 } 214 215 /* Restore vectors */ 216 PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x)); 217 PetscCall(DMDAVecRestoreArray(da, F, &f)); 218 PetscCall(DMRestoreLocalVector(da, &Xloc)); 219 PetscFunctionReturn(0); 220 } 221 222 /* --------------------------------------------------------------------- */ 223 /* 224 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 225 */ 226 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr) { 227 User user = (User)ptr; 228 DMDALocalInfo info; 229 PetscInt i; 230 DM da; 231 const Field *x, *xdot; 232 233 PetscFunctionBeginUser; 234 PetscCall(TSGetDM(ts, &da)); 235 PetscCall(DMDAGetLocalInfo(da, &info)); 236 237 /* Get pointers to vector data */ 238 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x)); 239 PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot)); 240 241 /* Compute function over the locally owned part of the grid */ 242 for (i = info.xs; i < info.xs + info.xm; i++) { 243 const PetscReal *k = user->k; 244 PetscScalar v[2][2]; 245 246 v[0][0] = a + k[0]; 247 v[0][1] = -k[1]; 248 v[1][0] = -k[0]; 249 v[1][1] = a + k[1]; 250 PetscCall(MatSetValuesBlocked(Jpre, 1, &i, 1, &i, &v[0][0], INSERT_VALUES)); 251 } 252 253 /* Restore vectors */ 254 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x)); 255 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot)); 256 257 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 258 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 259 if (J != Jpre) { 260 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 261 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 262 } 263 PetscFunctionReturn(0); 264 } 265 266 PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx) { 267 User user = (User)ctx; 268 DM da; 269 PetscInt i; 270 DMDALocalInfo info; 271 Field *x; 272 PetscReal hx; 273 274 PetscFunctionBeginUser; 275 PetscCall(TSGetDM(ts, &da)); 276 PetscCall(DMDAGetLocalInfo(da, &info)); 277 hx = 1.0 / (PetscReal)info.mx; 278 279 /* Get pointers to vector data */ 280 PetscCall(DMDAVecGetArray(da, X, &x)); 281 282 /* Compute function over the locally owned part of the grid */ 283 for (i = info.xs; i < info.xs + info.xm; i++) { 284 PetscReal r = (i + 1) * hx, ik = user->k[1] != 0.0 ? 1.0 / user->k[1] : 1.0; 285 x[i][0] = 1 + user->s[1] * r; 286 x[i][1] = user->k[0] * ik * x[i][0] + user->s[1] * ik; 287 } 288 PetscCall(DMDAVecRestoreArray(da, X, &x)); 289 PetscFunctionReturn(0); 290 } 291 292 /*TEST 293 294 test: 295 args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_arkimex_type 4 -ts_adapt_type none -ts_dt .005 -ts_max_time .1 296 requires: !single 297 298 test: 299 suffix: 2 300 args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_dt 1e-3 -ts_adapt_type none -ts_dt .005 -ts_max_time .1 301 nsize: 2 302 303 test: 304 suffix: 3 305 args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type ra34pw2 -ts_dt 5e-3 -ts_max_time .1 -ts_adapt_type none 306 nsize: 2 307 308 test: 309 suffix: 4 310 args: -ts_type eimex -da_grid_x 200 -ts_eimex_order_adapt true -ts_dt 0.001 -ts_monitor -ts_max_steps 100 311 filter: sed "s/ITS/TIME/g" 312 nsize: 2 313 314 TEST*/ 315