1 static const char help[] = "Time-dependent PDE in 1d. Simplified from ex15.c for illustrating how to solve DAEs. \n"; 2 /* 3 u_t = uxx 4 0 < x < 1; 5 At t=0: u(x) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5)) < .125 6 u(x) = 0.0 if r >= .125 7 8 Boundary conditions: 9 Dirichlet BC: 10 At x=0, x=1, u = 0.0 11 12 Neumann BC: 13 At x=0, x=1: du(x,t)/dx = 0 14 15 mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor 16 ./ex17 -da_grid_x 40 -monitor_solution 17 ./ex17 -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable 18 ./ex17 -jac_type fd_coloring -da_grid_x 500 -boundary 1 19 ./ex17 -da_grid_x 100 -ts_type gl -ts_adapt_type none -ts_max_steps 2 20 */ 21 22 #include <petscdm.h> 23 #include <petscdmda.h> 24 #include <petscts.h> 25 26 typedef enum { 27 JACOBIAN_ANALYTIC, 28 JACOBIAN_FD_COLORING, 29 JACOBIAN_FD_FULL 30 } JacobianType; 31 static const char *const JacobianTypes[] = {"analytic", "fd_coloring", "fd_full", "JacobianType", "fd_", 0}; 32 33 /* 34 User-defined data structures and routines 35 */ 36 typedef struct { 37 PetscReal c; 38 PetscInt boundary; /* Type of boundary condition */ 39 PetscBool viewJacobian; 40 } AppCtx; 41 42 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *); 43 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 44 static PetscErrorCode FormInitialSolution(TS, Vec, void *); 45 46 int main(int argc, char **argv) 47 { 48 TS ts; /* nonlinear solver */ 49 Vec u; /* solution, residual vectors */ 50 Mat J; /* Jacobian matrix */ 51 PetscInt nsteps; 52 PetscReal vmin, vmax, norm; 53 DM da; 54 PetscReal ftime, dt; 55 AppCtx user; /* user-defined work context */ 56 JacobianType jacType; 57 58 PetscFunctionBeginUser; 59 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 60 61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 Create distributed array (DMDA) to manage parallel grid and vectors 63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 64 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 1, 1, NULL, &da)); 65 PetscCall(DMSetFromOptions(da)); 66 PetscCall(DMSetUp(da)); 67 68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 Extract global vectors from DMDA; 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 PetscCall(DMCreateGlobalVector(da, &u)); 72 73 /* Initialize user application context */ 74 user.c = -30.0; 75 user.boundary = 0; /* 0: Dirichlet BC; 1: Neumann BC */ 76 user.viewJacobian = PETSC_FALSE; 77 78 PetscCall(PetscOptionsGetInt(NULL, NULL, "-boundary", &user.boundary, NULL)); 79 PetscCall(PetscOptionsHasName(NULL, NULL, "-viewJacobian", &user.viewJacobian)); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Create timestepping solver context 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 85 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 86 PetscCall(TSSetType(ts, TSTHETA)); 87 PetscCall(TSThetaSetTheta(ts, 1.0)); /* Make the Theta method behave like backward Euler */ 88 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user)); 89 90 PetscCall(DMSetMatType(da, MATAIJ)); 91 PetscCall(DMCreateMatrix(da, &J)); 92 jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */ 93 94 PetscCall(TSSetDM(ts, da)); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */ 95 96 ftime = 1.0; 97 PetscCall(TSSetMaxTime(ts, ftime)); 98 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Set initial conditions 102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 103 PetscCall(FormInitialSolution(ts, u, &user)); 104 PetscCall(TSSetSolution(ts, u)); 105 dt = .01; 106 PetscCall(TSSetTimeStep(ts, dt)); 107 108 /* Use slow fd Jacobian or fast fd Jacobian with colorings. 109 Note: this requirs snes which is not created until TSSetUp()/TSSetFromOptions() is called */ 110 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Options for Jacobian evaluation", NULL); 111 PetscCall(PetscOptionsEnum("-jac_type", "Type of Jacobian", "", JacobianTypes, (PetscEnum)jacType, (PetscEnum *)&jacType, 0)); 112 PetscOptionsEnd(); 113 if (jacType == JACOBIAN_ANALYTIC) { 114 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user)); 115 } else if (jacType == JACOBIAN_FD_COLORING) { 116 SNES snes; 117 PetscCall(TSGetSNES(ts, &snes)); 118 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, 0)); 119 } else if (jacType == JACOBIAN_FD_FULL) { 120 SNES snes; 121 PetscCall(TSGetSNES(ts, &snes)); 122 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, &user)); 123 } 124 125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 126 Set runtime options 127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 128 PetscCall(TSSetFromOptions(ts)); 129 130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 Integrate ODE system 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 PetscCall(TSSolve(ts, u)); 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Compute diagnostics of the solution 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 PetscCall(VecNorm(u, NORM_1, &norm)); 139 PetscCall(VecMax(u, NULL, &vmax)); 140 PetscCall(VecMin(u, NULL, &vmin)); 141 PetscCall(TSGetStepNumber(ts, &nsteps)); 142 PetscCall(TSGetTime(ts, &ftime)); 143 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "timestep %" PetscInt_FMT ": time %g, solution norm %g, max %g, min %g\n", nsteps, (double)ftime, (double)norm, (double)vmax, (double)vmin)); 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Free work space. 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 PetscCall(MatDestroy(&J)); 149 PetscCall(VecDestroy(&u)); 150 PetscCall(TSDestroy(&ts)); 151 PetscCall(DMDestroy(&da)); 152 PetscCall(PetscFinalize()); 153 return 0; 154 } 155 /* ------------------------------------------------------------------- */ 156 static PetscErrorCode FormIFunction(TS ts, PetscReal ftime, Vec U, Vec Udot, Vec F, void *ptr) 157 { 158 AppCtx *user = (AppCtx *)ptr; 159 DM da; 160 PetscInt i, Mx, xs, xm; 161 PetscReal hx, sx; 162 PetscScalar *u, *udot, *f; 163 Vec localU; 164 165 PetscFunctionBeginUser; 166 PetscCall(TSGetDM(ts, &da)); 167 PetscCall(DMGetLocalVector(da, &localU)); 168 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 169 170 hx = 1.0 / (PetscReal)(Mx - 1); 171 sx = 1.0 / (hx * hx); 172 173 /* 174 Scatter ghost points to local vector,using the 2-step process 175 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 176 By placing code between these two statements, computations can be 177 done while messages are in transition. 178 */ 179 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); 180 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); 181 182 /* Get pointers to vector data */ 183 PetscCall(DMDAVecGetArrayRead(da, localU, &u)); 184 PetscCall(DMDAVecGetArrayRead(da, Udot, &udot)); 185 PetscCall(DMDAVecGetArray(da, F, &f)); 186 187 /* Get local grid boundaries */ 188 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 189 190 /* Compute function over the locally owned part of the grid */ 191 for (i = xs; i < xs + xm; i++) { 192 if (user->boundary == 0) { /* Dirichlet BC */ 193 if (i == 0 || i == Mx - 1) f[i] = u[i]; /* F = U */ 194 else f[i] = udot[i] + (2. * u[i] - u[i - 1] - u[i + 1]) * sx; 195 } else { /* Neumann BC */ 196 if (i == 0) f[i] = u[0] - u[1]; 197 else if (i == Mx - 1) f[i] = u[i] - u[i - 1]; 198 else f[i] = udot[i] + (2. * u[i] - u[i - 1] - u[i + 1]) * sx; 199 } 200 } 201 202 /* Restore vectors */ 203 PetscCall(DMDAVecRestoreArrayRead(da, localU, &u)); 204 PetscCall(DMDAVecRestoreArrayRead(da, Udot, &udot)); 205 PetscCall(DMDAVecRestoreArray(da, F, &f)); 206 PetscCall(DMRestoreLocalVector(da, &localU)); 207 PetscFunctionReturn(0); 208 } 209 210 /* --------------------------------------------------------------------- */ 211 /* 212 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 213 */ 214 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *ctx) 215 { 216 PetscInt i, rstart, rend, Mx; 217 PetscReal hx, sx; 218 AppCtx *user = (AppCtx *)ctx; 219 DM da; 220 MatStencil col[3], row; 221 PetscInt nc; 222 PetscScalar vals[3]; 223 224 PetscFunctionBeginUser; 225 PetscCall(TSGetDM(ts, &da)); 226 PetscCall(MatGetOwnershipRange(Jpre, &rstart, &rend)); 227 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 228 hx = 1.0 / (PetscReal)(Mx - 1); 229 sx = 1.0 / (hx * hx); 230 for (i = rstart; i < rend; i++) { 231 nc = 0; 232 row.i = i; 233 if (user->boundary == 0 && (i == 0 || i == Mx - 1)) { 234 col[nc].i = i; 235 vals[nc++] = 1.0; 236 } else if (user->boundary > 0 && i == 0) { /* Left Neumann */ 237 col[nc].i = i; 238 vals[nc++] = 1.0; 239 col[nc].i = i + 1; 240 vals[nc++] = -1.0; 241 } else if (user->boundary > 0 && i == Mx - 1) { /* Right Neumann */ 242 col[nc].i = i - 1; 243 vals[nc++] = -1.0; 244 col[nc].i = i; 245 vals[nc++] = 1.0; 246 } else { /* Interior */ 247 col[nc].i = i - 1; 248 vals[nc++] = -1.0 * sx; 249 col[nc].i = i; 250 vals[nc++] = 2.0 * sx + a; 251 col[nc].i = i + 1; 252 vals[nc++] = -1.0 * sx; 253 } 254 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, vals, INSERT_VALUES)); 255 } 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 if (user->viewJacobian) { 264 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Jpre:\n")); 265 PetscCall(MatView(Jpre, PETSC_VIEWER_STDOUT_WORLD)); 266 } 267 PetscFunctionReturn(0); 268 } 269 270 /* ------------------------------------------------------------------- */ 271 PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ptr) 272 { 273 AppCtx *user = (AppCtx *)ptr; 274 PetscReal c = user->c; 275 DM da; 276 PetscInt i, xs, xm, Mx; 277 PetscScalar *u; 278 PetscReal hx, x, r; 279 280 PetscFunctionBeginUser; 281 PetscCall(TSGetDM(ts, &da)); 282 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 283 284 hx = 1.0 / (PetscReal)(Mx - 1); 285 286 /* Get pointers to vector data */ 287 PetscCall(DMDAVecGetArray(da, U, &u)); 288 289 /* Get local grid boundaries */ 290 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 291 292 /* Compute function over the locally owned part of the grid */ 293 for (i = xs; i < xs + xm; i++) { 294 x = i * hx; 295 r = PetscSqrtReal((x - .5) * (x - .5)); 296 if (r < .125) u[i] = PetscExpReal(c * r * r * r); 297 else u[i] = 0.0; 298 } 299 300 /* Restore vectors */ 301 PetscCall(DMDAVecRestoreArray(da, U, &u)); 302 PetscFunctionReturn(0); 303 } 304 305 /*TEST 306 307 test: 308 requires: !single 309 args: -da_grid_x 40 -ts_max_steps 2 -snes_monitor_short -ksp_monitor_short -ts_monitor 310 311 test: 312 suffix: 2 313 requires: !single 314 args: -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5 315 316 TEST*/ 317