https://mooseframework.inl.gov
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 
28 PetscErrorCode
29 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  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 
75 PetscErrorCode
76 PETScExternalSolverDestroy(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  */
87 PetscErrorCode
88 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  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 */
159 PetscErrorCode
160 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  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 */
295 PetscErrorCode
296 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  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 /* ------------------------------------------------------------------- */
409 PetscErrorCode
410 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  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 }
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
y
const std::vector< double > y
Definition: EquilibriumConstantFitTest.C:18
NS::F
static const std::string F
Definition: NS.h:165
NS::FV::converged
bool converged(const std::vector< std::pair< unsigned int, Real >> &residuals, const std::vector< Real > &abs_tolerances)
Based on the residuals, determine if the iterative process converged or not.
Definition: SegregatedSolverUtils.C:229
PETScExternalSolverCreate
PetscErrorCode PETScExternalSolverCreate(MPI_Comm comm, TS *ts)
Definition: PETScDiffusionFDM.C:29
x
const std::vector< double > x
Definition: EquilibriumConstantFitTest.C:17
f
Real f(Real x)
Test function for Brents method.
Definition: BrentsMethodTest.C:21
PETScExternalSolverDestroy
PetscErrorCode PETScExternalSolverDestroy(TS ts)
Definition: PETScDiffusionFDM.C:76
LibmeshPetscCallQ
LibmeshPetscCallQ(DMShellGetContext(dm, &ctx))
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
EM::j
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
PetscFunctionReturn
PetscFunctionReturn(LIBMESH_PETSC_SUCCESS)
Generated on Wed Aug 20 2025 08:33:58 for https://mooseframework.inl.gov by   doxygen 1.8.14