Hi,
we've been looking at using the built-in Schur complement capabilities in
PETSc, and encountered some issues which I think are bugs, though maybe I
just can't figure out how to use things properly.
Attached is a somewhat boiled down test program. It's based on
Crank-Nicolson TS, which uses SNES, which uses KSP, which has a PC that's
set to type fieldsplit. We'd like to have the KSP solver for the actual
Schur complement to use multigrid as PC, eventually.
One, kinda minor issue: when ComputeJacobian is called, we have to use
MatZeroEntries and recompute the Jacobian (of course, that's not really
needed for the simple problem we're solving right now because it's actually
linear). The reason for that is that, I think, the TS CN has shifted /
added a diagonal to the Jacobian matrix, as needed for CN. For this, I'm
just wondering if that's expected behavior (it's not obvious that one
wouldn't be passed the same matrix that was calculated the previous time).
Secondly, running the code with "-ts_view_pre" crashes, because the TS
isn't fully set up (in particular, the fieldsplit PC doesn't seem to have
properly created its sub-KSPs yet). I'm pretty sure it shouldn't crash
because of that. But maybe more importantly, the way things are set up
doesn't seem quite right. In particular, even after calling TSSetUp(), the
sub objects in the TS aren't properly set up yet -- that only seems to
happen in the first call to TSSolve() (or TSStep()). Trying to set up the
SNES inside the TS explicitly by calling SNESSetUp() doesn't work, either.
During TSStep, SNESSolve() is called with explicit vectors for b, x, it
appears only then TSSetUp() will actually work without errors.
I'm not sure when the fieldsplit PC should be set up to have its proper
sub-KSP, I suppose at least at TSSetUp() time, if not earlier. For us,
since we want to change the sub-KSPs, we need some way to do that before
the actual solve, but at that point, PCFieldSplitGetSubKSP() will fail --
this part of the code is #if 0'd out right now.
I believe that all of this works properly when using SNES directly, but we
can't get it to work inside the TS, so any help would be appreciated.
This is on the "maint" branch from petsc's git repo, from a couple of days
ago.
Thanks,
--Kai
--
Kai Germaschewski
Assistant Professor, Dept of Physics / Space Science Center
University of New Hampshire, Durham, NH 03824
office: Morse Hall 245E
phone: +1-603-862-2912
fax: +1-603-862-2771
static char help[] = "Schur.\n\n";
#include <petscdmda.h>
#include <petscts.h>
#include <stdio.h>
#undef __FUNCT__
#define __FUNCT__ __func__
#define CE CHKERRQ(ierr)
typedef struct {
PetscScalar u;
PetscScalar v;
} Field;
typedef struct {
Vec solution,rhs; /* global exact solution vector */
PetscInt m; /* total number of grid points */
PetscReal h; /* mesh width h = 2/(m-1) */
DM da;
} AppCtx;
// ----------------------------------------------------------------------
// SetInitialCondition
PetscErrorCode SetInitialCondition(AppCtx *appctx)
{
PetscInt xm,xs,i;
PetscReal h = appctx->h;
DM da = appctx->da;
Field *x;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = DMDAGetCorners(da, &xs,NULL, NULL, &xm,NULL, NULL); CE;
ierr = DMDAVecGetArray(da,appctx->solution, &x); CE;
for (i=xs; i<xs+xm; i++) {
x[i].u = PetscSinScalar(2*PETSC_PI*(i*h));
x[i].v = PetscSinScalar(2*PETSC_PI*(i*h));
}
ierr = DMDAVecRestoreArray(da,appctx->solution, &x); CE;
PetscFunctionReturn(0);
}
// ----------------------------------------------------------------------
// RHSFunction
PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx)
{
PetscErrorCode ierr;
Field *f,*xl;
AppCtx *appctx = (AppCtx*)ctx;
PetscInt i,xs,xm;
PetscReal h = appctx->h;
Vec XL;
DM da = appctx->da;
ierr = DMGetLocalVector(da,&XL); CE;
ierr = DMGlobalToLocalBegin(da, X, INSERT_VALUES, XL); CE;
ierr = DMGlobalToLocalEnd(da, X, INSERT_VALUES, XL); CE;
ierr = DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL); CE;
ierr = DMDAVecGetArray(da, F, &f); CE;
ierr = DMDAVecGetArray(da, XL, &xl); CE;
for (i=xs; i<xs+xm; i++) {
f[i].u = (xl[i+1].v - xl[i-1].v)/2/h;
f[i].v = (xl[i+1].u - xl[i-1].u)/2/h;
}
ierr = DMDAVecRestoreArray(da,F,&f); CE;
ierr = DMDAVecRestoreArray(da, XL, &xl); CE;
ierr = DMRestoreLocalVector(da,&XL); CE;
return(0);
}
// ----------------------------------------------------------------------
// ComputeJacobian
PetscErrorCode ComputeJacobian(TS ts,PetscReal t, Vec X, Mat *BB, Mat *A, MatStructure *flag, void *ctx)
{
AppCtx *appctx = (AppCtx*)ctx;
DM da = appctx->da;
PetscScalar h = appctx->h;
PetscErrorCode ierr;
PetscInt i,xm,xs;
PetscScalar v[2];
MatStencil row,col[2];
Mat B = *BB;
PetscScalar a = -1./(2*h);
ierr = MatZeroEntries(B); CE;
ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0); CE;
for (i=xs; i<xs+xm; i++) {
row.i = i; row.c = 0;
v[0] = a; col[0].i = i-1; col[0].c = 1;
v[1] = -a; col[1].i = i+1; col[1].c = 1;
ierr = MatSetValuesStencil(B,1,&row,2,col,v,INSERT_VALUES); CE;
row.c = 1;
v[0] = a; col[0].i = i-1; col[0].c = 0;
v[1] = -a; col[1].i = i+1; col[1].c = 0;
ierr = MatSetValuesStencil(B,1,&row,2,col,v,INSERT_VALUES); CE;
}
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CE;
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CE;
//*flag = SAME_NONZERO_PATTERN;
*flag = SAME_PRECONDITIONER;
PetscFunctionReturn(0);
}
// ----------------------------------------------------------------------
// DMDAVecOutput
PetscErrorCode DMDAVecOutput(DM dm, Vec xv, const char *name, int n_ghosts,int step)
{
PetscErrorCode ierr;
int ix,M;
const double Lx = 1;//2. * M_PI;
int xs, xm,ys;
char filename[strlen(name) + 20];
Field *x;
double xx;
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
sprintf(filename, "rpp%s_%06d_%06d.asc", name, step, rank);
FILE *f = fopen(filename, "w");
ierr = DMDAGetInfo(dm, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,NULL, NULL, NULL); CE;
ierr = DMDAGetCorners(dm, &xs, &ys, NULL, &xm, NULL, NULL); CE;
ierr = DMDAVecGetArray(dm, xv, &x); CE;
for (ix = xs; ix < xs + xm ; ix++) {
xx = Lx * ix / M;
fprintf(f, "%g %g %g\n", xx, x[ix].u, x[ix].v);
}
ierr = DMDAVecRestoreArray(dm, xv, &x); CE;
fclose(f);
PetscFunctionReturn(0);
}
// ----------------------------------------------------------------------
// Monitor
PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec b,void *ctx)
{
AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
PetscErrorCode ierr;
int out_every = 1;
if (step % out_every == 0) {
ierr = DMDAVecOutput(appctx->da, b, "x", 0, step); CE;
}
return 0;
}
// ----------------------------------------------------------------------
// main
int main(int argc,char **argv)
{
AppCtx appctx;
DM da;
TS ts;
PetscErrorCode ierr;
PetscInt step,m=16, time_steps_max = 20;
PetscReal dt = 0.1, time_total_max = 3.0;
ierr = PetscInitialize(&argc,&argv,(char*)0,help); CE;
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL); CE;
ierr = PetscOptionsGetInt(NULL,"-steps",&time_steps_max,NULL); CE;
appctx.m = m;
appctx.h = 1.0/appctx.m;
ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,
m,2,1,0,&da); CE;
appctx.da = da;
// state vector
ierr = DMGetGlobalVector(da,&appctx.solution); CE;
// set initial condition
ierr = SetInitialCondition(&appctx); CE;
ierr = TSCreate(PETSC_COMM_SELF,&ts); CE;
ierr = TSSetType(ts,TSCN); CE;
ierr = TSSetInitialTimeStep(ts, 0., dt); CE;
ierr = TSSetDuration(ts,time_steps_max,time_total_max); CE;
ierr = TSMonitorSet(ts,Monitor,&appctx,NULL); CE;
ierr = TSSetDM(ts,da); CE;
ierr = TSSetSolution(ts,appctx.solution); CE;
Mat J;
ierr = DMCreateMatrix(da,MATAIJ,&J); CE;
ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&appctx); CE;
ierr = TSSetRHSJacobian(ts,J,J,ComputeJacobian,&appctx); CE;
SNES snes;
ierr = TSGetSNES(ts, &snes); CE;
KSP ksp;
ierr = SNESGetKSP(snes, &ksp); CE;
PC pc;
ierr = KSPGetPC(ksp, &pc); CE;
PetscInt rs;
ierr = MatGetOwnershipRange(J,&rs,NULL); CE;
IS is[2];
ierr = ISCreateStride(PETSC_COMM_WORLD,m,rs,2,&is[0]); CE;
ierr = ISCreateStride(PETSC_COMM_WORLD,m,rs+1,2,&is[1]); CE;
ierr = PCSetType(pc,PCFIELDSPLIT); CE;
ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR); CE;
ierr = PCFieldSplitSetBlockSize(pc,2); CE;
ierr = PCFieldSplitSetIS(pc,NULL,is[0]); CE;
ierr = PCFieldSplitSetIS(pc,NULL,is[1]); CE;
ierr = ISDestroy(is); CE;
#if 0
int n;
KSP *subksp;
ierr = PCFieldSplitGetSubKSP(pc, &n, &subksp); CE;
PC subpc;
ierr = KSPGetPC(subksp[1], &subpc); CE;
ierr = PCSetType(subpc, PCMG); CE;
ierr = PetscFree(subksp); CE;
#endif
ierr = TSSetFromOptions(ts); CE;
ierr = TSSetUp(ts); CE;
ierr = TSSolve(ts,NULL); CE;
ierr = TSGetTimeStepNumber(ts,&step); CE;
ierr = TSDestroy(&ts); CE;
ierr = DMRestoreGlobalVector(da,&appctx.solution); CE;
ierr = MatDestroy(&J); CE;
ierr = DMDestroy(&da); CE;
ierr = PetscFinalize();
return 0;
}
