www.mooseframework.org
Functions
PETScDiffusionFDM.C File Reference

Go to the source code of this file.

Functions

PetscErrorCode PETScExternalSolverCreate (MPI_Comm comm, TS *ts)
 
PetscErrorCode PETScExternalSolverDestroy (TS ts)
 
PetscErrorCode externalPETScDiffusionFDMSolve (TS ts, Vec u0, Vec u, PetscReal dt, PetscReal time, PetscBool *converged)
 
PetscErrorCode FormIFunction (TS ts, PetscReal, Vec U, Vec Udot, Vec F, void *)
 
PetscErrorCode FormIJacobian (TS ts, PetscReal, Vec, Vec, PetscReal a, Mat J, Mat Jpre, void *)
 
PetscErrorCode FormInitialSolution (TS ts, Vec U, void *)
 

Function Documentation

◆ externalPETScDiffusionFDMSolve()

PetscErrorCode externalPETScDiffusionFDMSolve ( TS  ts,
Vec  u0,
Vec  u,
PetscReal  dt,
PetscReal  time,
PetscBool *  converged 
)

Definition at line 105 of file PETScDiffusionFDM.C.

Referenced by ExternalPETScProblem::externalSolve().

107 {
108  PetscErrorCode ierr;
109  TSConvergedReason reason;
110 #if !PETSC_VERSION_LESS_THAN(3, 8, 0)
111  PetscInt current_step;
112 #endif
113  DM da;
114 
115  PetscFunctionBeginUser;
116  PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
117  PetscValidType(ts, 1);
118  PetscValidHeaderSpecific(u0, VEC_CLASSID, 2);
119  PetscValidType(u0, 2);
120  PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
121  PetscValidType(u, 3);
122 #if PETSC_VERSION_LESS_THAN(3, 19, 3)
123  PetscValidPointer(converged, 6);
124 #else
125  PetscAssertPointer(converged, 6);
126 #endif
127 
128  ierr = TSGetDM(ts, &da);
129  CHKERRQ(ierr);
130 
131 #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
132  PetscOptionsSetValue(NULL, "-ts_monitor", NULL);
133  PetscOptionsSetValue(NULL, "-snes_monitor", NULL);
134  PetscOptionsSetValue(NULL, "-ksp_monitor", NULL);
135 #else
136  PetscOptionsSetValue("-ts_monitor", NULL);
137  PetscOptionsSetValue("-snes_monitor", NULL);
138  PetscOptionsSetValue("-ksp_monitor", NULL);
139 #endif
140 
141  /*ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr);*/
142  ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
143  CHKERRQ(ierr);
144 
145  ierr = VecCopy(u0, u);
146  CHKERRQ(ierr);
147 
148  ierr = TSSetSolution(ts, u);
149  CHKERRQ(ierr);
150  ierr = TSSetTimeStep(ts, dt);
151  CHKERRQ(ierr);
152  ierr = TSSetTime(ts, time - dt);
153  CHKERRQ(ierr);
154 #if !PETSC_VERSION_LESS_THAN(3, 8, 0)
155  ierr = TSGetStepNumber(ts, &current_step);
156  CHKERRQ(ierr);
157  ierr = TSSetMaxSteps(ts, current_step + 1);
158  CHKERRQ(ierr);
159 #else
160  SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Require PETSc-3.8.x or higher ");
161 #endif
162  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163  Sets various TS parameters from user options
164  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165  ierr = TSSetFromOptions(ts);
166  CHKERRQ(ierr);
167  /*
168  * Let the passed-in dt take the priority
169  */
170  ierr = TSSetTimeStep(ts, dt);
171  CHKERRQ(ierr);
172  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173  Solve nonlinear system
174  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175  ierr = TSSolve(ts, u);
176  CHKERRQ(ierr);
177 
178  ierr = TSGetConvergedReason(ts, &reason);
179  CHKERRQ(ierr);
180  *converged = reason > 0 ? PETSC_TRUE : PETSC_FALSE;
181 
182  PetscFunctionReturn(0);
183 }
ierr
CHKERRQ(ierr)

◆ FormIFunction()

PetscErrorCode FormIFunction ( TS  ts,
PetscReal  ,
Vec  U,
Vec  Udot,
Vec  F,
void  
)

Definition at line 190 of file PETScDiffusionFDM.C.

Referenced by PETScExternalSolverCreate().

191 {
192  PetscErrorCode ierr;
193  DM da;
194  PetscInt i, j, Mx, My, xs, ys, xm, ym;
195  PetscReal hx, hy, sx, sy;
196  PetscScalar u, uxx, uyy, **uarray, **f, **udot;
197  Vec localU;
198  MPI_Comm comm;
199 
200  PetscFunctionBeginUser;
201  ierr = PetscObjectGetComm((PetscObject)ts, &comm);
202  ierr = TSGetDM(ts, &da);
203  CHKERRQ(ierr);
204  ierr = DMGetLocalVector(da, &localU);
205  CHKERRQ(ierr);
206  ierr = DMDAGetInfo(da,
207  PETSC_IGNORE,
208  &Mx,
209  &My,
210  PETSC_IGNORE,
211  PETSC_IGNORE,
212  PETSC_IGNORE,
213  PETSC_IGNORE,
214  PETSC_IGNORE,
215  PETSC_IGNORE,
216  PETSC_IGNORE,
217  PETSC_IGNORE,
218  PETSC_IGNORE,
219  PETSC_IGNORE);
220  CHKERRQ(ierr);
221 
222  hx = 1.0 / (PetscReal)(Mx - 1);
223  sx = 1.0 / (hx * hx);
224  hy = 1.0 / (PetscReal)(My - 1);
225  sy = 1.0 / (hy * hy);
226 
227  /*
228  Scatter ghost points to local vector,using the 2-step process
229  DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
230  By placing code between these two statements, computations can be
231  done while messages are in transition.
232  */
233  ierr = DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
234  CHKERRQ(ierr);
235  ierr = DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
236  CHKERRQ(ierr);
237 
238  /* Get pointers to vector data */
239  ierr = DMDAVecGetArrayRead(da, localU, &uarray);
240  CHKERRQ(ierr);
241  ierr = DMDAVecGetArray(da, F, &f);
242  CHKERRQ(ierr);
243  ierr = DMDAVecGetArray(da, Udot, &udot);
244  CHKERRQ(ierr);
245 
246  /* Get local grid boundaries */
247  ierr = DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
248  CHKERRQ(ierr);
249 
250  /* Compute function over the locally owned part of the grid */
251  for (j = ys; j < ys + ym; j++)
252  {
253  for (i = xs; i < xs + xm; i++)
254  {
255  /* Boundary conditions */
256  if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1)
257  {
258  if (PETSC_TRUE)
259  { /* Drichlet BC */
260  f[j][i] = uarray[j][i]; /* F = U */
261  }
262  else
263  { /* Neumann BC */
264  if (i == 0 && j == 0)
265  { /* SW corner */
266  f[j][i] = uarray[j][i] - uarray[j + 1][i + 1];
267  }
268  else if (i == Mx - 1 && j == 0)
269  { /* SE corner */
270  f[j][i] = uarray[j][i] - uarray[j + 1][i - 1];
271  }
272  else if (i == 0 && j == My - 1)
273  { /* NW corner */
274  f[j][i] = uarray[j][i] - uarray[j - 1][i + 1];
275  }
276  else if (i == Mx - 1 && j == My - 1)
277  { /* NE corner */
278  f[j][i] = uarray[j][i] - uarray[j - 1][i - 1];
279  }
280  else if (i == 0)
281  { /* Left */
282  f[j][i] = uarray[j][i] - uarray[j][i + 1];
283  }
284  else if (i == Mx - 1)
285  { /* Right */
286  f[j][i] = uarray[j][i] - uarray[j][i - 1];
287  }
288  else if (j == 0)
289  { /* Bottom */
290  f[j][i] = uarray[j][i] - uarray[j + 1][i];
291  }
292  else if (j == My - 1)
293  { /* Top */
294  f[j][i] = uarray[j][i] - uarray[j - 1][i];
295  }
296  }
297  }
298  else
299  { /* Interior */
300  u = uarray[j][i];
301  /* 5-point stencil */
302  uxx = (-2.0 * u + uarray[j][i - 1] + uarray[j][i + 1]);
303  uyy = (-2.0 * u + uarray[j - 1][i] + uarray[j + 1][i]);
304  if (PETSC_FALSE)
305  {
306  /* 9-point stencil: assume hx=hy */
307  uxx = 2.0 * uxx / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] +
308  uarray[j + 1][i - 1] + uarray[j + 1][i + 1]) -
309  2.0 * u) /
310  6.0;
311  uyy = 2.0 * uyy / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] +
312  uarray[j + 1][i - 1] + uarray[j + 1][i + 1]) -
313  2.0 * u) /
314  6.0;
315  }
316  f[j][i] = udot[j][i] - (uxx * sx + uyy * sy);
317  }
318  }
319  }
320 
321  /* Restore vectors */
322  ierr = DMDAVecRestoreArrayRead(da, localU, &uarray);
323  CHKERRQ(ierr);
324  ierr = DMDAVecRestoreArray(da, F, &f);
325  CHKERRQ(ierr);
326  ierr = DMDAVecRestoreArray(da, Udot, &udot);
327  CHKERRQ(ierr);
328  ierr = DMRestoreLocalVector(da, &localU);
329  CHKERRQ(ierr);
330  ierr = PetscLogFlops(11.0 * ym * xm);
331  CHKERRQ(ierr);
332  PetscFunctionReturn(0);
333 }
ierr
static const std::string F
Definition: NS.h:150
Real f(Real x)
Test function for Brents method.
CHKERRQ(ierr)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ FormIJacobian()

PetscErrorCode FormIJacobian ( TS  ts,
PetscReal  ,
Vec  ,
Vec  ,
PetscReal  a,
Mat  J,
Mat  Jpre,
void  
)

Definition at line 341 of file PETScDiffusionFDM.C.

Referenced by PETScExternalSolverCreate().

343 {
344  PetscErrorCode ierr;
345  PetscInt i, j, Mx, My, xs, ys, xm, ym, nc;
346  DM da;
347  MatStencil col[5], row;
348  PetscScalar vals[5], hx, hy, sx, sy;
349 
350  PetscFunctionBeginUser;
351  ierr = TSGetDM(ts, &da);
352  CHKERRQ(ierr);
353  ierr = DMDAGetInfo(da,
354  PETSC_IGNORE,
355  &Mx,
356  &My,
357  PETSC_IGNORE,
358  PETSC_IGNORE,
359  PETSC_IGNORE,
360  PETSC_IGNORE,
361  PETSC_IGNORE,
362  PETSC_IGNORE,
363  PETSC_IGNORE,
364  PETSC_IGNORE,
365  PETSC_IGNORE,
366  PETSC_IGNORE);
367  CHKERRQ(ierr);
368  ierr = DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
369  CHKERRQ(ierr);
370 
371  hx = 1.0 / (PetscReal)(Mx - 1);
372  sx = 1.0 / (hx * hx);
373  hy = 1.0 / (PetscReal)(My - 1);
374  sy = 1.0 / (hy * hy);
375 
376  for (j = ys; j < ys + ym; j++)
377  {
378  for (i = xs; i < xs + xm; i++)
379  {
380  nc = 0;
381  row.j = j;
382  row.i = i;
383  if (PETSC_TRUE && (i == 0 || i == Mx - 1 || j == 0 || j == My - 1))
384  {
385  col[nc].j = j;
386  col[nc].i = i;
387  vals[nc++] = 1.0;
388  }
389  else if (PETSC_FALSE && i == 0)
390  { /* Left Neumann */
391  col[nc].j = j;
392  col[nc].i = i;
393  vals[nc++] = 1.0;
394  col[nc].j = j;
395  col[nc].i = i + 1;
396  vals[nc++] = -1.0;
397  }
398  else if (PETSC_FALSE && i == Mx - 1)
399  { /* Right Neumann */
400  col[nc].j = j;
401  col[nc].i = i;
402  vals[nc++] = 1.0;
403  col[nc].j = j;
404  col[nc].i = i - 1;
405  vals[nc++] = -1.0;
406  }
407  else if (PETSC_FALSE && j == 0)
408  { /* Bottom Neumann */
409  col[nc].j = j;
410  col[nc].i = i;
411  vals[nc++] = 1.0;
412  col[nc].j = j + 1;
413  col[nc].i = i;
414  vals[nc++] = -1.0;
415  }
416  else if (PETSC_FALSE && j == My - 1)
417  { /* Top Neumann */
418  col[nc].j = j;
419  col[nc].i = i;
420  vals[nc++] = 1.0;
421  col[nc].j = j - 1;
422  col[nc].i = i;
423  vals[nc++] = -1.0;
424  }
425  else
426  { /* Interior */
427  col[nc].j = j - 1;
428  col[nc].i = i;
429  vals[nc++] = -sy;
430  col[nc].j = j;
431  col[nc].i = i - 1;
432  vals[nc++] = -sx;
433  col[nc].j = j;
434  col[nc].i = i;
435  vals[nc++] = 2.0 * (sx + sy) + a;
436  col[nc].j = j;
437  col[nc].i = i + 1;
438  vals[nc++] = -sx;
439  col[nc].j = j + 1;
440  col[nc].i = i;
441  vals[nc++] = -sy;
442  }
443  ierr = MatSetValuesStencil(Jpre, 1, &row, nc, col, vals, INSERT_VALUES);
444  CHKERRQ(ierr);
445  }
446  }
447  ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
448  CHKERRQ(ierr);
449  ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
450  CHKERRQ(ierr);
451  if (J != Jpre)
452  {
453  ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
454  CHKERRQ(ierr);
455  ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
456  CHKERRQ(ierr);
457  }
458 
459  PetscFunctionReturn(0);
460 }
ierr
CHKERRQ(ierr)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ FormInitialSolution()

PetscErrorCode FormInitialSolution ( TS  ts,
Vec  U,
void  
)

Definition at line 464 of file PETScDiffusionFDM.C.

Referenced by ExternalPETScProblem::ExternalPETScProblem().

465 {
466  DM da;
467  PetscReal c = -30.0;
468  PetscErrorCode ierr;
469  PetscInt i, j, xs, ys, xm, ym, Mx, My;
470  PetscScalar ** u;
471  PetscReal hx, hy, x, y, r;
472 
473  PetscFunctionBeginUser;
474  ierr = TSGetDM(ts, &da);
475  CHKERRQ(ierr);
476  ierr = DMDAGetInfo(da,
477  PETSC_IGNORE,
478  &Mx,
479  &My,
480  PETSC_IGNORE,
481  PETSC_IGNORE,
482  PETSC_IGNORE,
483  PETSC_IGNORE,
484  PETSC_IGNORE,
485  PETSC_IGNORE,
486  PETSC_IGNORE,
487  PETSC_IGNORE,
488  PETSC_IGNORE,
489  PETSC_IGNORE);
490  CHKERRQ(ierr);
491 
492  hx = 1.0 / (PetscReal)(Mx - 1);
493  hy = 1.0 / (PetscReal)(My - 1);
494 
495  /* Get pointers to vector data */
496  ierr = DMDAVecGetArray(da, U, &u);
497  CHKERRQ(ierr);
498 
499  /* Get local grid boundaries */
500  ierr = DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
501  CHKERRQ(ierr);
502 
503  /* Compute function over the locally owned part of the grid */
504  for (j = ys; j < ys + ym; j++)
505  {
506  y = j * hy;
507  for (i = xs; i < xs + xm; i++)
508  {
509  x = i * hx;
510  r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
511  if (r < .125)
512  u[j][i] = PetscExpReal(c * r * r * r);
513  else
514  u[j][i] = 0.0;
515  }
516  }
517 
518  /* Restore vectors */
519  ierr = DMDAVecRestoreArray(da, U, &u);
520  CHKERRQ(ierr);
521  PetscFunctionReturn(0);
522 }
const std::vector< double > y
ierr
const std::vector< double > x
CHKERRQ(ierr)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ PETScExternalSolverCreate()

PetscErrorCode PETScExternalSolverCreate ( MPI_Comm  comm,
TS *  ts 
)

Definition at line 29 of file PETScDiffusionFDM.C.

Referenced by ExternalPetscSolverApp::getPetscTS().

30 {
31  DM da;
32 
33  PetscErrorCode ierr;
34 
35  PetscFunctionBeginUser;
36  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37  Create distributed array (DMDA) to manage parallel grid and vectors
38  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39  ierr = DMDACreate2d(comm,
40  DM_BOUNDARY_NONE,
41  DM_BOUNDARY_NONE,
42  DMDA_STENCIL_STAR,
43  11,
44  11,
45  PETSC_DECIDE,
46  PETSC_DECIDE,
47  1,
48  1,
49  NULL,
50  NULL,
51  &da);
52  CHKERRQ(ierr);
53  ierr = DMSetFromOptions(da);
54  CHKERRQ(ierr);
55  ierr = DMSetUp(da);
56  CHKERRQ(ierr);
57 
58  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59  Create timestepping solver context
60  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61  ierr = TSCreate(comm, ts);
62  CHKERRQ(ierr);
63  ierr = TSSetProblemType(*ts, TS_NONLINEAR);
64  CHKERRQ(ierr);
65  ierr = TSSetType(*ts, TSBEULER);
66  CHKERRQ(ierr);
67  ierr = TSSetDM(*ts, da);
68  CHKERRQ(ierr);
69  ierr = DMDestroy(&da);
70  CHKERRQ(ierr);
71  ierr = TSSetIFunction(*ts, NULL, FormIFunction, nullptr);
72  CHKERRQ(ierr);
73  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74  Set Jacobian evaluation routine
75  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76  ierr = TSSetIJacobian(*ts, NULL, NULL, FormIJacobian, NULL);
77  CHKERRQ(ierr);
78  /*
79  * Make it consistent with the original solver
80  * This can be changed during runtime via -ts_dt
81  */
82  ierr = TSSetTimeStep(*ts, 0.2);
83  CHKERRQ(ierr);
84  ierr = TSSetFromOptions(*ts);
85  CHKERRQ(ierr);
86  PetscFunctionReturn(0);
87 }
PetscErrorCode FormIJacobian(TS ts, PetscReal, Vec, Vec, PetscReal a, Mat J, Mat Jpre, void *)
PetscErrorCode FormIFunction(TS ts, PetscReal, Vec U, Vec Udot, Vec F, void *)
ierr
CHKERRQ(ierr)

◆ PETScExternalSolverDestroy()

PetscErrorCode PETScExternalSolverDestroy ( TS  ts)

Definition at line 90 of file PETScDiffusionFDM.C.

Referenced by ExternalPetscSolverApp::~ExternalPetscSolverApp().

91 {
92  PetscErrorCode ierr;
93 
94  PetscFunctionBeginUser;
95  ierr = TSDestroy(&ts);
96  CHKERRQ(ierr);
97  PetscFunctionReturn(0);
98 }
ierr
CHKERRQ(ierr)