https://mooseframework.inl.gov
Loading...
Searching...
No Matches
PETScDiffusionFDM.C
Go to the documentation of this file.
1//* This file is part of the MOOSE framework
2//* https://mooseframework.inl.gov
3//*
4//* All rights reserved, see COPYRIGHT for full restrictions
5//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6//*
7//* Licensed under LGPL 2.1, please see LICENSE for details
8//* https://www.gnu.org/licenses/lgpl-2.1.html
9
10#include <PETScDiffusionFDM.h>
11
12/*
13 u_t = uxx + uyy
14 0 < x < 1, 0 < y < 1;
15 At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125
16 u(x,y) = 0.0 if r >= .125
17
18
19 Boundary conditions:
20 Drichlet BC:
21 At x=0, x=1, y=0, y=1: u = 0.0
22
23 Neumann BC:
24 At x=0, x=1: du(x,y,t)/dx = 0
25 At y=0, y=1: du(x,y,t)/dy = 0
26*/
27
28PetscErrorCode
29PETScExternalSolverCreate(MPI_Comm comm, TS * ts)
30{
31 DM da;
32
33 PetscFunctionBeginUser;
34 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35 Create distributed array (DMDA) to manage parallel grid and vectors
36 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37 LibmeshPetscCallQ(DMDACreate2d(comm,
38 DM_BOUNDARY_NONE,
39 DM_BOUNDARY_NONE,
40 DMDA_STENCIL_STAR,
41 11,
42 11,
43 PETSC_DECIDE,
44 PETSC_DECIDE,
45 1,
46 1,
47 NULL,
48 NULL,
49 &da));
50 LibmeshPetscCallQ(DMSetFromOptions(da));
51 LibmeshPetscCallQ(DMSetUp(da));
52
53 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54 Create timestepping solver context
55 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56 LibmeshPetscCallQ(TSCreate(comm, ts));
57 LibmeshPetscCallQ(TSSetProblemType(*ts, TS_NONLINEAR));
58 LibmeshPetscCallQ(TSSetType(*ts, TSBEULER));
59 LibmeshPetscCallQ(TSSetDM(*ts, da));
60 LibmeshPetscCallQ(DMDestroy(&da));
61 LibmeshPetscCallQ(TSSetIFunction(*ts, NULL, FormIFunction, nullptr));
62 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 Set Jacobian evaluation routine
64 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 LibmeshPetscCallQ(TSSetIJacobian(*ts, NULL, NULL, FormIJacobian, NULL));
66 /*
67 * Make it consistent with the original solver
68 * This can be changed during runtime via -ts_dt
69 */
70 LibmeshPetscCallQ(TSSetTimeStep(*ts, 0.2));
71 LibmeshPetscCallQ(TSSetFromOptions(*ts));
72 PetscFunctionReturn(PETSC_SUCCESS);
73}
74
75PetscErrorCode
76PETScExternalSolverDestroy(TS ts)
77{
78 PetscFunctionBeginUser;
79 LibmeshPetscCallQ(TSDestroy(&ts));
80 PetscFunctionReturn(PETSC_SUCCESS);
81}
82
83/*
84 * This is a modified version of PETSc/src/ts/examples/tutorials/ex15.c
85 * to demonstrate how MOOSE interact with an external solver package
86 */
87PetscErrorCode
88externalPETScDiffusionFDMSolve(
89 TS ts, Vec u0, Vec u, PetscReal dt, PetscReal time, PetscBool * converged)
90{
91 TSConvergedReason reason;
92#if !PETSC_VERSION_LESS_THAN(3, 8, 0)
93 PetscInt current_step;
94#endif
95 DM da;
96
97 PetscFunctionBeginUser;
98 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
99 PetscValidType(ts, 1);
100 PetscValidHeaderSpecific(u0, VEC_CLASSID, 2);
101 PetscValidType(u0, 2);
102 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
103 PetscValidType(u, 3);
104#if PETSC_VERSION_LESS_THAN(3, 19, 3)
105 PetscValidPointer(converged, 6);
106#else
107 PetscAssertPointer(converged, 6);
108#endif
109
110 LibmeshPetscCallQ(TSGetDM(ts, &da));
111
112#if !PETSC_VERSION_LESS_THAN(3, 7, 0)
113 LibmeshPetscCallQ(PetscOptionsSetValue(NULL, "-ts_monitor", NULL));
114 LibmeshPetscCallQ(PetscOptionsSetValue(NULL, "-snes_monitor", NULL));
115 LibmeshPetscCallQ(PetscOptionsSetValue(NULL, "-ksp_monitor", NULL));
116#else
117 LibmeshPetscCallQ(PetscOptionsSetValue("-ts_monitor", NULL));
118 LibmeshPetscCallQ(PetscOptionsSetValue("-snes_monitor", NULL));
119 LibmeshPetscCallQ(PetscOptionsSetValue("-ksp_monitor", NULL));
120#endif
121
122 /*ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr);*/
123 LibmeshPetscCallQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
124
125 LibmeshPetscCallQ(VecCopy(u0, u));
126
127 LibmeshPetscCallQ(TSSetSolution(ts, u));
128 LibmeshPetscCallQ(TSSetTimeStep(ts, dt));
129 LibmeshPetscCallQ(TSSetTime(ts, time - dt));
130#if !PETSC_VERSION_LESS_THAN(3, 8, 0)
131 LibmeshPetscCallQ(TSGetStepNumber(ts, &current_step));
132 LibmeshPetscCallQ(TSSetMaxSteps(ts, current_step + 1));
133#else
134 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Require PETSc-3.8.x or higher ");
135#endif
136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137 Sets various TS parameters from user options
138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139 LibmeshPetscCallQ(TSSetFromOptions(ts));
140 /*
141 * Let the passed-in dt take the priority
142 */
143 LibmeshPetscCallQ(TSSetTimeStep(ts, dt));
144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145 Solve nonlinear system
146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147 LibmeshPetscCallQ(TSSolve(ts, u));
148
149 LibmeshPetscCallQ(TSGetConvergedReason(ts, &reason));
150 *converged = reason > 0 ? PETSC_TRUE : PETSC_FALSE;
151
152 PetscFunctionReturn(PETSC_SUCCESS);
153}
154
155/* --------------------------------------------------------------------- */
156/*
157 FormIFunction = Udot - RHSFunction
158*/
159PetscErrorCode
160FormIFunction(TS ts, PetscReal /*t*/, Vec U, Vec Udot, Vec F, void * /*ctx*/)
161{
162 DM da;
163 PetscInt i, j, Mx, My, xs, ys, xm, ym;
164 PetscReal hx, hy, sx, sy;
165 PetscScalar u, uxx, uyy, **uarray, **f, **udot;
166 Vec localU;
167 MPI_Comm comm;
168
169 PetscFunctionBeginUser;
170 LibmeshPetscCallQ(PetscObjectGetComm((PetscObject)ts, &comm));
171 LibmeshPetscCallQ(TSGetDM(ts, &da));
172 LibmeshPetscCallQ(DMGetLocalVector(da, &localU));
173 LibmeshPetscCallQ(DMDAGetInfo(da,
174 PETSC_IGNORE,
175 &Mx,
176 &My,
177 PETSC_IGNORE,
178 PETSC_IGNORE,
179 PETSC_IGNORE,
180 PETSC_IGNORE,
181 PETSC_IGNORE,
182 PETSC_IGNORE,
183 PETSC_IGNORE,
184 PETSC_IGNORE,
185 PETSC_IGNORE,
186 PETSC_IGNORE));
187
188 hx = 1.0 / (PetscReal)(Mx - 1);
189 sx = 1.0 / (hx * hx);
190 hy = 1.0 / (PetscReal)(My - 1);
191 sy = 1.0 / (hy * hy);
192
193 /*
194 Scatter ghost points to local vector,using the 2-step process
195 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
196 By placing code between these two statements, computations can be
197 done while messages are in transition.
198 */
199 LibmeshPetscCallQ(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
200 LibmeshPetscCallQ(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
201
202 /* Get pointers to vector data */
203 LibmeshPetscCallQ(DMDAVecGetArrayRead(da, localU, &uarray));
204 LibmeshPetscCallQ(DMDAVecGetArray(da, F, &f));
205 LibmeshPetscCallQ(DMDAVecGetArray(da, Udot, &udot));
206
207 /* Get local grid boundaries */
208 LibmeshPetscCallQ(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
209
210 /* Compute function over the locally owned part of the grid */
211 for (j = ys; j < ys + ym; j++)
212 {
213 for (i = xs; i < xs + xm; i++)
214 {
215 /* Boundary conditions */
216 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1)
217 {
218 if (PETSC_TRUE)
219 { /* Drichlet BC */
220 f[j][i] = uarray[j][i]; /* F = U */
221 }
222 else
223 { /* Neumann BC */
224 if (i == 0 && j == 0)
225 { /* SW corner */
226 f[j][i] = uarray[j][i] - uarray[j + 1][i + 1];
227 }
228 else if (i == Mx - 1 && j == 0)
229 { /* SE corner */
230 f[j][i] = uarray[j][i] - uarray[j + 1][i - 1];
231 }
232 else if (i == 0 && j == My - 1)
233 { /* NW corner */
234 f[j][i] = uarray[j][i] - uarray[j - 1][i + 1];
235 }
236 else if (i == Mx - 1 && j == My - 1)
237 { /* NE corner */
238 f[j][i] = uarray[j][i] - uarray[j - 1][i - 1];
239 }
240 else if (i == 0)
241 { /* Left */
242 f[j][i] = uarray[j][i] - uarray[j][i + 1];
243 }
244 else if (i == Mx - 1)
245 { /* Right */
246 f[j][i] = uarray[j][i] - uarray[j][i - 1];
247 }
248 else if (j == 0)
249 { /* Bottom */
250 f[j][i] = uarray[j][i] - uarray[j + 1][i];
251 }
252 else if (j == My - 1)
253 { /* Top */
254 f[j][i] = uarray[j][i] - uarray[j - 1][i];
255 }
256 }
257 }
258 else
259 { /* Interior */
260 u = uarray[j][i];
261 /* 5-point stencil */
262 uxx = (-2.0 * u + uarray[j][i - 1] + uarray[j][i + 1]);
263 uyy = (-2.0 * u + uarray[j - 1][i] + uarray[j + 1][i]);
264 if (PETSC_FALSE)
265 {
266 /* 9-point stencil: assume hx=hy */
267 uxx = 2.0 * uxx / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] +
268 uarray[j + 1][i - 1] + uarray[j + 1][i + 1]) -
269 2.0 * u) /
270 6.0;
271 uyy = 2.0 * uyy / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] +
272 uarray[j + 1][i - 1] + uarray[j + 1][i + 1]) -
273 2.0 * u) /
274 6.0;
275 }
276 f[j][i] = udot[j][i] - (uxx * sx + uyy * sy);
277 }
278 }
279 }
280
281 /* Restore vectors */
282 LibmeshPetscCallQ(DMDAVecRestoreArrayRead(da, localU, &uarray));
283 LibmeshPetscCallQ(DMDAVecRestoreArray(da, F, &f));
284 LibmeshPetscCallQ(DMDAVecRestoreArray(da, Udot, &udot));
285 LibmeshPetscCallQ(DMRestoreLocalVector(da, &localU));
286 LibmeshPetscCallQ(PetscLogFlops(11.0 * ym * xm));
287 PetscFunctionReturn(PETSC_SUCCESS);
288}
289
290/* --------------------------------------------------------------------- */
291/*
292 FormIJacobian() - Compute IJacobian = dF/dU + a dF/dUdot
293 This routine is not used with option '-use_coloring'
294*/
295PetscErrorCode
296FormIJacobian(
297 TS ts, PetscReal /*t*/, Vec /*U*/, Vec /*Udot*/, PetscReal a, Mat J, Mat Jpre, void * /*ctx*/)
298{
299 PetscInt i, j, Mx, My, xs, ys, xm, ym, nc;
300 DM da;
301 MatStencil col[5], row;
302 PetscScalar vals[5], hx, hy, sx, sy;
303
304 PetscFunctionBeginUser;
305 LibmeshPetscCallQ(TSGetDM(ts, &da));
306 LibmeshPetscCallQ(DMDAGetInfo(da,
307 PETSC_IGNORE,
308 &Mx,
309 &My,
310 PETSC_IGNORE,
311 PETSC_IGNORE,
312 PETSC_IGNORE,
313 PETSC_IGNORE,
314 PETSC_IGNORE,
315 PETSC_IGNORE,
316 PETSC_IGNORE,
317 PETSC_IGNORE,
318 PETSC_IGNORE,
319 PETSC_IGNORE));
320 LibmeshPetscCallQ(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
321
322 hx = 1.0 / (PetscReal)(Mx - 1);
323 sx = 1.0 / (hx * hx);
324 hy = 1.0 / (PetscReal)(My - 1);
325 sy = 1.0 / (hy * hy);
326
327 for (j = ys; j < ys + ym; j++)
328 {
329 for (i = xs; i < xs + xm; i++)
330 {
331 nc = 0;
332 row.j = j;
333 row.i = i;
334 if (PETSC_TRUE && (i == 0 || i == Mx - 1 || j == 0 || j == My - 1))
335 {
336 col[nc].j = j;
337 col[nc].i = i;
338 vals[nc++] = 1.0;
339 }
340 else if (PETSC_FALSE && i == 0)
341 { /* Left Neumann */
342 col[nc].j = j;
343 col[nc].i = i;
344 vals[nc++] = 1.0;
345 col[nc].j = j;
346 col[nc].i = i + 1;
347 vals[nc++] = -1.0;
348 }
349 else if (PETSC_FALSE && i == Mx - 1)
350 { /* Right Neumann */
351 col[nc].j = j;
352 col[nc].i = i;
353 vals[nc++] = 1.0;
354 col[nc].j = j;
355 col[nc].i = i - 1;
356 vals[nc++] = -1.0;
357 }
358 else if (PETSC_FALSE && j == 0)
359 { /* Bottom Neumann */
360 col[nc].j = j;
361 col[nc].i = i;
362 vals[nc++] = 1.0;
363 col[nc].j = j + 1;
364 col[nc].i = i;
365 vals[nc++] = -1.0;
366 }
367 else if (PETSC_FALSE && j == My - 1)
368 { /* Top Neumann */
369 col[nc].j = j;
370 col[nc].i = i;
371 vals[nc++] = 1.0;
372 col[nc].j = j - 1;
373 col[nc].i = i;
374 vals[nc++] = -1.0;
375 }
376 else
377 { /* Interior */
378 col[nc].j = j - 1;
379 col[nc].i = i;
380 vals[nc++] = -sy;
381 col[nc].j = j;
382 col[nc].i = i - 1;
383 vals[nc++] = -sx;
384 col[nc].j = j;
385 col[nc].i = i;
386 vals[nc++] = 2.0 * (sx + sy) + a;
387 col[nc].j = j;
388 col[nc].i = i + 1;
389 vals[nc++] = -sx;
390 col[nc].j = j + 1;
391 col[nc].i = i;
392 vals[nc++] = -sy;
393 }
394 LibmeshPetscCallQ(MatSetValuesStencil(Jpre, 1, &row, nc, col, vals, INSERT_VALUES));
395 }
396 }
397 LibmeshPetscCallQ(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
398 LibmeshPetscCallQ(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
399 if (J != Jpre)
400 {
401 LibmeshPetscCallQ(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
402 LibmeshPetscCallQ(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
403 }
404
405 PetscFunctionReturn(PETSC_SUCCESS);
406}
407
408/* ------------------------------------------------------------------- */
409PetscErrorCode
410FormInitialSolution(TS ts, Vec U, void * /*ptr*/)
411{
412 DM da;
413 PetscReal c = -30.0;
414 PetscInt i, j, xs, ys, xm, ym, Mx, My;
415 PetscScalar ** u;
416 PetscReal hx, hy, x, y, r;
417
418 PetscFunctionBeginUser;
419 LibmeshPetscCallQ(TSGetDM(ts, &da));
420 LibmeshPetscCallQ(DMDAGetInfo(da,
421 PETSC_IGNORE,
422 &Mx,
423 &My,
424 PETSC_IGNORE,
425 PETSC_IGNORE,
426 PETSC_IGNORE,
427 PETSC_IGNORE,
428 PETSC_IGNORE,
429 PETSC_IGNORE,
430 PETSC_IGNORE,
431 PETSC_IGNORE,
432 PETSC_IGNORE,
433 PETSC_IGNORE));
434
435 hx = 1.0 / (PetscReal)(Mx - 1);
436 hy = 1.0 / (PetscReal)(My - 1);
437
438 /* Get pointers to vector data */
439 LibmeshPetscCallQ(DMDAVecGetArray(da, U, &u));
440
441 /* Get local grid boundaries */
442 LibmeshPetscCallQ(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
443
444 /* Compute function over the locally owned part of the grid */
445 for (j = ys; j < ys + ym; j++)
446 {
447 y = j * hy;
448 for (i = xs; i < xs + xm; i++)
449 {
450 x = i * hx;
451 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
452 if (r < .125)
453 u[j][i] = PetscExpReal(c * r * r * r);
454 else
455 u[j][i] = 0.0;
456 }
457 }
458
459 /* Restore vectors */
460 LibmeshPetscCallQ(DMDAVecRestoreArray(da, U, &u));
461 PetscFunctionReturn(PETSC_SUCCESS);
462}
f
Real f(Real x)
Test function for Brents method.
Definition BrentsMethodTest.C:21
y
const std::vector< double > y
Definition EquilibriumConstantFitTest.C:18
x
const std::vector< double > x
Definition EquilibriumConstantFitTest.C:17
PETScExternalSolverDestroy
PetscErrorCode PETScExternalSolverDestroy(TS ts)
Definition PETScDiffusionFDM.C:76
FormInitialSolution
PetscErrorCode FormInitialSolution(TS ts, Vec U, void *)
Definition PETScDiffusionFDM.C:410
externalPETScDiffusionFDMSolve
PetscErrorCode externalPETScDiffusionFDMSolve(TS ts, Vec u0, Vec u, PetscReal dt, PetscReal time, PetscBool *converged)
Definition PETScDiffusionFDM.C:88
PETScExternalSolverCreate
PetscErrorCode PETScExternalSolverCreate(MPI_Comm comm, TS *ts)
Definition PETScDiffusionFDM.C:29
FormIJacobian
PetscErrorCode FormIJacobian(TS ts, PetscReal, Vec, Vec, PetscReal a, Mat J, Mat Jpre, void *)
Definition PETScDiffusionFDM.C:296
FormIFunction
PetscErrorCode FormIFunction(TS ts, PetscReal, Vec U, Vec Udot, Vec F, void *)
Definition PETScDiffusionFDM.C:160
Generated on Fri Aug 21 2026 13:40:59 for https://mooseframework.inl.gov by doxygen 1.9.8