Line data Source code
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 :
28 : PetscErrorCode
29 370 : PETScExternalSolverCreate(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 370 : 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 370 : LibmeshPetscCallQ(DMSetFromOptions(da));
51 370 : LibmeshPetscCallQ(DMSetUp(da));
52 :
53 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54 : Create timestepping solver context
55 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56 370 : LibmeshPetscCallQ(TSCreate(comm, ts));
57 370 : LibmeshPetscCallQ(TSSetProblemType(*ts, TS_NONLINEAR));
58 370 : LibmeshPetscCallQ(TSSetType(*ts, TSBEULER));
59 370 : LibmeshPetscCallQ(TSSetDM(*ts, da));
60 370 : LibmeshPetscCallQ(DMDestroy(&da));
61 370 : LibmeshPetscCallQ(TSSetIFunction(*ts, NULL, FormIFunction, nullptr));
62 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 : Set Jacobian evaluation routine
64 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 370 : 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 370 : LibmeshPetscCallQ(TSSetTimeStep(*ts, 0.2));
71 370 : LibmeshPetscCallQ(TSSetFromOptions(*ts));
72 370 : PetscFunctionReturn(PETSC_SUCCESS);
73 : }
74 :
75 : PetscErrorCode
76 749 : PETScExternalSolverDestroy(TS ts)
77 : {
78 : PetscFunctionBeginUser;
79 749 : LibmeshPetscCallQ(TSDestroy(&ts));
80 749 : 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 : */
87 : PetscErrorCode
88 5925 : externalPETScDiffusionFDMSolve(
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 5925 : LibmeshPetscCallQ(TSGetDM(ts, &da));
111 :
112 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
113 5925 : LibmeshPetscCallQ(PetscOptionsSetValue(NULL, "-ts_monitor", NULL));
114 5925 : LibmeshPetscCallQ(PetscOptionsSetValue(NULL, "-snes_monitor", NULL));
115 5925 : 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 5925 : LibmeshPetscCallQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
124 :
125 5925 : LibmeshPetscCallQ(VecCopy(u0, u));
126 :
127 5925 : LibmeshPetscCallQ(TSSetSolution(ts, u));
128 5925 : LibmeshPetscCallQ(TSSetTimeStep(ts, dt));
129 5925 : LibmeshPetscCallQ(TSSetTime(ts, time - dt));
130 : #if !PETSC_VERSION_LESS_THAN(3, 8, 0)
131 5925 : LibmeshPetscCallQ(TSGetStepNumber(ts, ¤t_step));
132 5925 : 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 5925 : LibmeshPetscCallQ(TSSetFromOptions(ts));
140 : /*
141 : * Let the passed-in dt take the priority
142 : */
143 5925 : LibmeshPetscCallQ(TSSetTimeStep(ts, dt));
144 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145 : Solve nonlinear system
146 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147 5925 : LibmeshPetscCallQ(TSSolve(ts, u));
148 :
149 5925 : LibmeshPetscCallQ(TSGetConvergedReason(ts, &reason));
150 5925 : *converged = reason > 0 ? PETSC_TRUE : PETSC_FALSE;
151 :
152 5925 : PetscFunctionReturn(PETSC_SUCCESS);
153 : }
154 :
155 : /* --------------------------------------------------------------------- */
156 : /*
157 : FormIFunction = Udot - RHSFunction
158 : */
159 : PetscErrorCode
160 17775 : FormIFunction(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 17775 : LibmeshPetscCallQ(PetscObjectGetComm((PetscObject)ts, &comm));
171 17775 : LibmeshPetscCallQ(TSGetDM(ts, &da));
172 17775 : LibmeshPetscCallQ(DMGetLocalVector(da, &localU));
173 17775 : 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 17775 : hx = 1.0 / (PetscReal)(Mx - 1);
189 17775 : sx = 1.0 / (hx * hx);
190 17775 : hy = 1.0 / (PetscReal)(My - 1);
191 17775 : 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 17775 : LibmeshPetscCallQ(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
200 17775 : LibmeshPetscCallQ(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
201 :
202 : /* Get pointers to vector data */
203 17775 : LibmeshPetscCallQ(DMDAVecGetArrayRead(da, localU, &uarray));
204 17775 : LibmeshPetscCallQ(DMDAVecGetArray(da, F, &f));
205 17775 : LibmeshPetscCallQ(DMDAVecGetArray(da, Udot, &udot));
206 :
207 : /* Get local grid boundaries */
208 17775 : LibmeshPetscCallQ(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
209 :
210 : /* Compute function over the locally owned part of the grid */
211 124596 : for (j = ys; j < ys + ym; j++)
212 : {
213 1149720 : for (i = xs; i < xs + xm; i++)
214 : {
215 : /* Boundary conditions */
216 1042899 : if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1)
217 : {
218 344760 : if (PETSC_TRUE)
219 : { /* Drichlet BC */
220 344760 : 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 698139 : u = uarray[j][i];
261 : /* 5-point stencil */
262 698139 : uxx = (-2.0 * u + uarray[j][i - 1] + uarray[j][i + 1]);
263 698139 : 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 698139 : f[j][i] = udot[j][i] - (uxx * sx + uyy * sy);
277 : }
278 : }
279 : }
280 :
281 : /* Restore vectors */
282 17775 : LibmeshPetscCallQ(DMDAVecRestoreArrayRead(da, localU, &uarray));
283 17775 : LibmeshPetscCallQ(DMDAVecRestoreArray(da, F, &f));
284 17775 : LibmeshPetscCallQ(DMDAVecRestoreArray(da, Udot, &udot));
285 17775 : LibmeshPetscCallQ(DMRestoreLocalVector(da, &localU));
286 17775 : LibmeshPetscCallQ(PetscLogFlops(11.0 * ym * xm));
287 17775 : 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 : */
295 : PetscErrorCode
296 11850 : FormIJacobian(
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 11850 : LibmeshPetscCallQ(TSGetDM(ts, &da));
306 11850 : 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 11850 : LibmeshPetscCallQ(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
321 :
322 11850 : hx = 1.0 / (PetscReal)(Mx - 1);
323 11850 : sx = 1.0 / (hx * hx);
324 11850 : hy = 1.0 / (PetscReal)(My - 1);
325 11850 : sy = 1.0 / (hy * hy);
326 :
327 83064 : for (j = ys; j < ys + ym; j++)
328 : {
329 766480 : for (i = xs; i < xs + xm; i++)
330 : {
331 : nc = 0;
332 695266 : row.j = j;
333 695266 : row.i = i;
334 695266 : if (PETSC_TRUE && (i == 0 || i == Mx - 1 || j == 0 || j == My - 1))
335 : {
336 229840 : col[nc].j = j;
337 229840 : col[nc].i = i;
338 229840 : 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 465426 : col[nc].j = j - 1;
379 465426 : col[nc].i = i;
380 465426 : vals[nc++] = -sy;
381 465426 : col[nc].j = j;
382 465426 : col[nc].i = i - 1;
383 465426 : vals[nc++] = -sx;
384 465426 : col[nc].j = j;
385 465426 : col[nc].i = i;
386 465426 : vals[nc++] = 2.0 * (sx + sy) + a;
387 465426 : col[nc].j = j;
388 465426 : col[nc].i = i + 1;
389 465426 : vals[nc++] = -sx;
390 465426 : col[nc].j = j + 1;
391 465426 : col[nc].i = i;
392 465426 : vals[nc++] = -sy;
393 : }
394 695266 : LibmeshPetscCallQ(MatSetValuesStencil(Jpre, 1, &row, nc, col, vals, INSERT_VALUES));
395 : }
396 : }
397 11850 : LibmeshPetscCallQ(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
398 11850 : LibmeshPetscCallQ(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
399 11850 : if (J != Jpre)
400 : {
401 0 : LibmeshPetscCallQ(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
402 0 : LibmeshPetscCallQ(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
403 : }
404 :
405 11850 : PetscFunctionReturn(PETSC_SUCCESS);
406 : }
407 :
408 : /* ------------------------------------------------------------------- */
409 : PetscErrorCode
410 370 : FormInitialSolution(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 370 : LibmeshPetscCallQ(TSGetDM(ts, &da));
420 370 : 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 370 : hx = 1.0 / (PetscReal)(Mx - 1);
436 370 : hy = 1.0 / (PetscReal)(My - 1);
437 :
438 : /* Get pointers to vector data */
439 370 : LibmeshPetscCallQ(DMDAVecGetArray(da, U, &u));
440 :
441 : /* Get local grid boundaries */
442 370 : LibmeshPetscCallQ(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
443 :
444 : /* Compute function over the locally owned part of the grid */
445 2482 : for (j = ys; j < ys + ym; j++)
446 : {
447 2112 : y = j * hy;
448 22440 : for (i = xs; i < xs + xm; i++)
449 : {
450 20328 : x = i * hx;
451 20328 : r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
452 20328 : if (r < .125)
453 840 : u[j][i] = PetscExpReal(c * r * r * r);
454 : else
455 19488 : u[j][i] = 0.0;
456 : }
457 : }
458 :
459 : /* Restore vectors */
460 370 : LibmeshPetscCallQ(DMDAVecRestoreArray(da, U, &u));
461 370 : PetscFunctionReturn(PETSC_SUCCESS);
462 : }
|
