Lines Matching +full:- +full:j
1 static char help[] = "Time-dependent PDE in 2d. Modified from ex13.c for illustrating how to solve …
5 At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125
16 mpiexec -n 2 ./ex15 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
17 ./ex15 -da_grid_x 40 -da_grid_y 40 -draw_pause .1 -boundary 1 -ts_monitor_draw_solution
18 ./ex15 -da_grid_x 40 -da_grid_y 40 -draw_pause .1 -boundary 1 -Jtype 2 -nstencilpts 9
27 User-defined data structures and routines
47 Mat J, Jmf = NULL; /* Jacobian matrices */ in main() local
50 AppCtx user; /* user-defined work context */ in main()
62 user.c = -30.0; in main()
66 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nstencilpts", &user.nstencilpts, NULL)); in main()
67 PetscCall(PetscOptionsGetInt(NULL, NULL, "-boundary", &user.boundary, NULL)); in main()
68 PetscCall(PetscOptionsHasName(NULL, NULL, "-viewJacobian", &user.viewJacobian)); in main()
70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
90 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
111 PetscCall(DMCreateMatrix(da, &J)); in main()
113 PetscCall(PetscOptionsGetInt(NULL, NULL, "-Jtype", &Jtype, NULL)); in main()
116 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user)); in main()
117 …} else { /* use finite difference Jacobian J as preconditioner and '-snes_mf_operator' for Mat*vec… in main()
120 if (Jtype == 1) { /* slow finite difference J; */ in main()
121 PetscCall(SNESSetJacobian(snes, Jmf, J, SNESComputeJacobianDefault, NULL)); in main()
122 } else if (Jtype == 2) { /* Use coloring to compute finite difference J efficiently */ in main()
123 PetscCall(SNESSetJacobian(snes, Jmf, J, SNESComputeJacobianDefaultColor, 0)); in main()
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
132 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
134 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
140 PetscCall(MatDestroy(&J)); in main()
151 /* --------------------------------------------------------------------- */
153 FormIFunction = Udot - RHSFunction
158 DM da = user->da; in FormIFunction()
159 PetscInt i, j, Mx, My, xs, ys, xm, ym; in FormIFunction() local
168 hx = 1.0 / (PetscReal)(Mx - 1); in FormIFunction()
170 hy = 1.0 / (PetscReal)(My - 1); in FormIFunction()
172 …PetscCheck(user->nstencilpts != 9 || hx == hy, PETSC_COMM_WORLD, PETSC_ERR_SUP, "hx must equal hy … in FormIFunction()
175 Scatter ghost points to local vector,using the 2-step process in FormIFunction()
192 for (j = ys; j < ys + ym; j++) { in FormIFunction()
195 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) { in FormIFunction()
196 if (user->boundary == 0) { /* Drichlet BC */ in FormIFunction()
197 f[j][i] = uarray[j][i]; /* F = U */ in FormIFunction()
199 if (i == 0 && j == 0) { /* SW corner */ in FormIFunction()
200 f[j][i] = uarray[j][i] - uarray[j + 1][i + 1]; in FormIFunction()
201 } else if (i == Mx - 1 && j == 0) { /* SE corner */ in FormIFunction()
202 f[j][i] = uarray[j][i] - uarray[j + 1][i - 1]; in FormIFunction()
203 } else if (i == 0 && j == My - 1) { /* NW corner */ in FormIFunction()
204 f[j][i] = uarray[j][i] - uarray[j - 1][i + 1]; in FormIFunction()
205 } else if (i == Mx - 1 && j == My - 1) { /* NE corner */ in FormIFunction()
206 f[j][i] = uarray[j][i] - uarray[j - 1][i - 1]; in FormIFunction()
208 f[j][i] = uarray[j][i] - uarray[j][i + 1]; in FormIFunction()
209 } else if (i == Mx - 1) { /* Right */ in FormIFunction()
210 f[j][i] = uarray[j][i] - uarray[j][i - 1]; in FormIFunction()
211 } else if (j == 0) { /* Bottom */ in FormIFunction()
212 f[j][i] = uarray[j][i] - uarray[j + 1][i]; in FormIFunction()
213 } else if (j == My - 1) { /* Top */ in FormIFunction()
214 f[j][i] = uarray[j][i] - uarray[j - 1][i]; in FormIFunction()
218 u = uarray[j][i]; in FormIFunction()
219 /* 5-point stencil */ in FormIFunction()
220 uxx = (-2.0 * u + uarray[j][i - 1] + uarray[j][i + 1]); in FormIFunction()
221 uyy = (-2.0 * u + uarray[j - 1][i] + uarray[j + 1][i]); in FormIFunction()
222 if (user->nstencilpts == 9) { in FormIFunction()
223 /* 9-point stencil: assume hx=hy */ in FormIFunction()
224 … 2.0 * uxx / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] + uarray[j + 1][i - 1] + ua… in FormIFunction()
225 … 2.0 * uyy / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] + uarray[j + 1][i - 1] + ua… in FormIFunction()
227 f[j][i] = udot[j][i] - (uxx * sx + uyy * sy); in FormIFunction()
241 /* --------------------------------------------------------------------- */
243 FormIJacobian() - Compute IJacobian = dF/dU + a dF/dUdot
244 This routine is not used with option '-use_coloring'
246 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, Pet… in FormIJacobian() argument
248 PetscInt i, j, Mx, My, xs, ys, xm, ym, nc; in FormIJacobian() local
250 DM da = user->da; in FormIJacobian()
258 hx = 1.0 / (PetscReal)(Mx - 1); in FormIJacobian()
260 hy = 1.0 / (PetscReal)(My - 1); in FormIJacobian()
263 for (j = ys; j < ys + ym; j++) { in FormIJacobian()
266 row.j = j; in FormIJacobian()
268 if (user->boundary == 0 && (i == 0 || i == Mx - 1 || j == 0 || j == My - 1)) { in FormIJacobian()
269 col[nc].j = j; in FormIJacobian()
273 } else if (user->boundary > 0 && i == 0) { /* Left Neumann */ in FormIJacobian()
274 col[nc].j = j; in FormIJacobian()
277 col[nc].j = j; in FormIJacobian()
279 vals[nc++] = -1.0; in FormIJacobian()
280 } else if (user->boundary > 0 && i == Mx - 1) { /* Right Neumann */ in FormIJacobian()
281 col[nc].j = j; in FormIJacobian()
284 col[nc].j = j; in FormIJacobian()
285 col[nc].i = i - 1; in FormIJacobian()
286 vals[nc++] = -1.0; in FormIJacobian()
287 } else if (user->boundary > 0 && j == 0) { /* Bottom Neumann */ in FormIJacobian()
288 col[nc].j = j; in FormIJacobian()
291 col[nc].j = j + 1; in FormIJacobian()
293 vals[nc++] = -1.0; in FormIJacobian()
294 } else if (user->boundary > 0 && j == My - 1) { /* Top Neumann */ in FormIJacobian()
295 col[nc].j = j; in FormIJacobian()
298 col[nc].j = j - 1; in FormIJacobian()
300 vals[nc++] = -1.0; in FormIJacobian()
302 col[nc].j = j - 1; in FormIJacobian()
304 vals[nc++] = -sy; in FormIJacobian()
305 col[nc].j = j; in FormIJacobian()
306 col[nc].i = i - 1; in FormIJacobian()
307 vals[nc++] = -sx; in FormIJacobian()
308 col[nc].j = j; in FormIJacobian()
311 col[nc].j = j; in FormIJacobian()
313 vals[nc++] = -sx; in FormIJacobian()
314 col[nc].j = j + 1; in FormIJacobian()
316 vals[nc++] = -sy; in FormIJacobian()
323 if (J != Jpre) { in FormIJacobian()
324 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); in FormIJacobian()
325 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); in FormIJacobian()
328 if (user->viewJacobian) { in FormIJacobian()
335 /* ------------------------------------------------------------------- */
339 DM da = user->da; in FormInitialSolution()
340 PetscReal c = user->c; in FormInitialSolution()
341 PetscInt i, j, xs, ys, xm, ym, Mx, My; in FormInitialSolution() local
348 hx = 1.0 / (PetscReal)(Mx - 1); in FormInitialSolution()
349 hy = 1.0 / (PetscReal)(My - 1); in FormInitialSolution()
358 for (j = ys; j < ys + ym; j++) { in FormInitialSolution()
359 y = j * hy; in FormInitialSolution()
362 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5)); in FormInitialSolution()
363 if (r < .125) u[j][i] = PetscExpReal(c * r * r * r); in FormInitialSolution()
364 else u[j][i] = 0.0; in FormInitialSolution()
376 args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -ts_monitor
380 args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -Jtype 2 -ts_monitor
385 args: -da_grid_x 20 -da_grid_y 20 -boundary 1 -ts_max_steps 10 -ts_monitor
391 args: -da_grid_x 20 -da_grid_y 20 -boundary 1 -ts_max_steps 10 -ts_monitor
396 args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -Jtype 1 -ts_monitor