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