static char help[] = "Solve a toy 1D problem on a staggered grid.\n\
Accepts command line options -a, -b, and -c\n\
and approximately solves\n\
u''(x) = c, u(0) = 1, u(1) = b\n\n";
/*
To demonstrate the basic functionality of DMStag, solves a second-order ODE,
u''(x) = f(x), 0 < x < 1
u(0) = a
u(1) = b
in mixed form, that is by introducing an auxiliary variable p
p'(x) = f(x), p - u'(x) = 0, 0 < x < 1
u(0) = a
u(1) = b
For f == c, the solution is
u(x) = a + (b - a - (c/2)) x + (c/2) x^2
p(x) = (b - a - (c/2)) + c x
To find an approximate solution, discretize by storing values of p in
elements and values of u on their boundaries, and using first-order finite
differences.
This should in fact produce a (nodal) solution with no discretization error,
so differences from the reference solution will be restricted to those induced
by floating point operations (in particular, the finite differences) and the
accuracy of the linear solve.
Parameters for the main grid can be controlled with command line options, e.g.
-stag_grid_x 10
In particular to notice in this example are the two methods of indexing. The
first is analogous to the use of MatStencil with DMDA, and the second is
analogous to the use of DMDAVecGetArrayDOF().
The first, recommended for ease of use, is based on naming an element in the
global grid, a location in its support, and a component. For example,
one might consider element e, the left side (a vertex in 1d), and the first
component (index 0). This is accomplished by populating a DMStagStencil struct,
e.g.
DMStagStencil stencil;
stencil.i = i
stencil.loc = DMSTAG_LEFT;
stencil.c = 0
Note that below, for convenenience, we #define an alias LEFT for DMSTAG_LEFT.
The second, which ensures maximum efficiency, makes use of the underlying
block structure of local DMStag-derived vectors, and requires the user to
obtain the correct local offset for the degrees of freedom they would like to
use. This is made simple with the helper function DMStagGetLocationSlot().
Note that the linear system being solved is indefinite, so is not entirely
trivial to invert. The default solver here (GMRES/Jacobi) is a poor choice,
made to avoid depending on an external package. To solve a larger system,
the usual method for a 1-d problem such as this is to choose a sophisticated
direct solver, e.g. configure --download-suitesparse and run
$PETSC_DIR/$PETSC_ARCH/bin/mpiexec -n 3 ./stag_ex2 -stag_grid_x 100 -pc_type lu -pc_factor_mat_solver_package umfpack
You can also impose a periodic boundary condition, in which case -b and -c are
ignored; b = a and c = 0.0 are used, giving a constant u == a , p == 0.
-stag_boundary_type_x periodic
*/
#include
#include
#include
/* Shorter, more convenient names for DMStagStencilLocation entries */
#define LEFT DMSTAG_LEFT
#define RIGHT DMSTAG_RIGHT
#define ELEMENT DMSTAG_ELEMENT
int main(int argc,char **argv)
{
PetscErrorCode ierr;
DM dmSol,dmForcing;
DM dmCoordSol;
Vec sol,solRef,solRefLocal,f,fLocal,rhs,coordSolLocal;
Mat A;
PetscScalar a,b,c,h;
KSP ksp;
PC pc;
PetscInt start,n,e,nExtra;
PetscInt iu,ip,ixu,ixp;
PetscBool isLastRank,isFirstRank;
PetscScalar **arrSol,**arrCoordSol;
DMBoundaryType boundary;
const PetscReal domainSize = 1.0;
/* Initialize PETSc */
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
/* Create 1D DMStag for the solution, and set up. Note that you can supply many
command line options (see the man page for DMStagCreate1d)
*/
ierr = DMStagCreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,1,1,DMSTAG_STENCIL_BOX,1,NULL,&dmSol);CHKERRQ(ierr);
ierr = DMSetFromOptions(dmSol);CHKERRQ(ierr);
ierr = DMSetUp(dmSol);CHKERRQ(ierr);
/* Create uniform coordinates. Note that in higher-dimensional examples,
coordinates are created differently.*/
ierr = DMStagSetUniformCoordinatesExplicit(dmSol,0.0,domainSize,0.0,0.0,0.0,0.0);CHKERRQ(ierr);
/* Determine boundary type */
ierr = DMStagGetBoundaryTypes(dmSol,&boundary,NULL,NULL);CHKERRQ(ierr);
/* Process command line options (depends on DMStag setup) */
a = 1.0; b = 2.0; c = 1.0;
ierr = PetscOptionsGetScalar(NULL,NULL,"-a",&a,NULL);CHKERRQ(ierr);
if (boundary == DM_BOUNDARY_PERIODIC) {
b = a;
c = 0.0;
} else {
ierr = PetscOptionsGetScalar(NULL,NULL,"-b",&b,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetScalar(NULL,NULL,"-c",&c,NULL);CHKERRQ(ierr);
}
/* Compute reference solution on the grid, using direct array access */
ierr = DMCreateGlobalVector(dmSol,&solRef);CHKERRQ(ierr);
ierr = DMGetLocalVector(dmSol,&solRefLocal);CHKERRQ(ierr);
ierr = DMStagVecGetArray(dmSol,solRefLocal,&arrSol);CHKERRQ(ierr);
ierr = DMGetCoordinateDM(dmSol,&dmCoordSol);CHKERRQ(ierr);
ierr = DMGetCoordinatesLocal(dmSol,&coordSolLocal);CHKERRQ(ierr);
ierr = DMStagVecGetArrayRead(dmCoordSol,coordSolLocal,&arrCoordSol);CHKERRQ(ierr);
ierr = DMStagGetCorners(dmSol,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);CHKERRQ(ierr);
/* Get the correct entries for each of our variables in local element-wise storage */
ierr = DMStagGetLocationSlot(dmSol,LEFT, 0,&iu);CHKERRQ(ierr);
ierr = DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);CHKERRQ(ierr);
ierr = DMStagGetLocationSlot(dmCoordSol,LEFT, 0,&ixu);CHKERRQ(ierr);
ierr = DMStagGetLocationSlot(dmCoordSol,ELEMENT,0,&ixp);CHKERRQ(ierr);
for (e=start; e