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