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