Here,
run this e.g. with
mpirun -np 4 ./test1 -draw_pause -1 -pc_type mg -pc_mg_galerkin
In Debug mode it crashes, in standard mode it fails. Using Petsc current.
Run without mg and works!
Hope it's not me being stupid and doing something terribly wrong!
Best,
Filippo
On Monday 20 October 2014 07:29:04 Barry Smith wrote:
> Can you send us the code: to [email protected] or
> [email protected] ? Or something that reproduces the problem?
>
> Barry
>
> > On Oct 20, 2014, at 3:31 AM, Filippo Leonardi
> > <[email protected]> wrote:
> >
> > Hi,
> >
> > I have a very specific problem that I cannot figure out with PCMG and
> > multiple solves. I got a linear system that I solve many times, with same
> > matrix but different RHS. I can successfully solve the system with
> > standard techniques, such as default solver or LU or PCGAMG. Even MG
> > works if I destroy the ksp each time or I recompute the matrices at each
> > time. But when I try and go to MG and not recomputing the matrices each
> > time the solver fails. Any idea?
> >
> > Here some detail: the setup:
> > ierr = KSPSetDMActive(ksp, PETSC_FALSE); CHKERRQ(ierr);
> > ierr = KSPSetDM(ksp, da); CHKERRQ(ierr);
> > ierr = KSPSetComputeOperators(ksp, ComputeMatrix, ctx); CHKERRQ(ierr);
> > ierr = KSPSetComputeRHS(ksp, ComputeRHS, ctx); CHKERRQ(ierr);
> > ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
> >
> > Then I solve as usual, for a large number of time steps:
> > ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
> > ierr = KSPGetSolution(ksp, &phi);
> >
> > The solver converges and does that in a reasonable number of iterations:
> > Linear solve converged due to CONVERGED_RTOL iterations 7
> > And ksp_view and ksp_monitor do not show any weird stuff.
> >
> > - Weirdly enough using any solver (for instance cg+bjacobi or gamg)
> > everything works (So the matrix and RHS are working fine).
> > - But the problem persists with Galerkin matrices (-pc_mg_galerkin) (So is
> > not a ComputeMatrix problem).
> > - If I do:
> > ierr = KSPSetComputeOperators(ksp, ComputeMatrix, this); CHKERRQ(ierr);
> > between each solve or destroy the ksp entirely each time the solution is
> > also perfect (So is not a boundary scaling or other stuff problem).
> > - If I run MG with only 2 levels (so just coarse) I also get a fine
> > result.
> >
> > Setup RHS is called each time as expected, setup matrix is called just
> > once, also as expected.
> >
> > The only thing I can think is that MG does not update some value that
> > actually needs to be recomputed.
> >
> > Any idea?
> >
> > The solution is not that different from:
> > http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/
> > ex29.c.html
> >
> > Best,
> > Filippo
static char help[] = "Parallel vector layout.\n\n";
#include <petscsys.h>
#include <petscvec.h>
#include <petscdmda.h>
#include <petscksp.h>
DM da_sc;
KSP ksp;
PetscErrorCode ierr;
struct UserContext {
Vec rhs;
};
int N = 10;
#undef __FUNCT__
#define __FUNCT__ "ComputeRHS"
PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx) {
UserContext* user = (UserContext*) ctx;
DM da;
PetscInt mx, my;
PetscScalar Hx, Hy;
ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
Hx = 1. / (PetscReal) (mx);
Hy = 1. / (PetscReal) (my);
ierr = VecScale(user->rhs, Hx * Hy); CHKERRQ(ierr);
ierr = VecCopy(user->rhs, b); CHKERRQ(ierr);
MatNullSpace nullspace;
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace); CHKERRQ(ierr);
ierr = MatNullSpaceRemove(nullspace, b); CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&nullspace); CHKERRQ(ierr);
ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
return 0;
}
#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx) {
PetscInt mx, my, xm, ym, xs, ys;
PetscScalar Hx, Hy;
DM da;
MatStencil row;
PetscScalar HydHx, HxdHy;
PetscScalar v_c;
PetscInt ncols;
ierr = KSPGetDM(ksp, &da);CHKERRQ(ierr);
ierr = DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); CHKERRQ(ierr);
Hx = 1. / (PetscReal) (mx);
Hy = 1. / (PetscReal) (my);
ierr = MatSetOption(jac, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); CHKERRQ(ierr);
// ierr = MatSetOption(jac, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(ierr);
HxdHy = Hx/Hy;
HydHx = Hy/Hx;
ierr = DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0); CHKERRQ(ierr);
PetscScalar v[5];
int k = 0;
for (int j=ys; j<ys+ym; ++j) {
for (int i=xs; i<xs+xm; ++i) {
ncols = 0;
MatStencil col[5] = {{0}};
v_c = 0;
row.k = k; row.j = j; row.i = i; row.c = 0;
col[ncols].k=k; col[ncols].j = j; col[ncols].i = i-1; v[ncols++] = HydHx; v_c -= HydHx;
col[ncols].k=k; col[ncols].j = j; col[ncols].i = i+1; v[ncols++] = HydHx; v_c -= HydHx;
col[ncols].k=k; col[ncols].j = j-1; col[ncols].i = i; v[ncols++] = HxdHy; v_c -= HxdHy;
col[ncols].k=k; col[ncols].j = j+1; col[ncols].i = i; v[ncols++] = HxdHy; v_c -= HxdHy;
col[ncols].k=k; col[ncols].j=j; col[ncols].i = i; v[ncols++] = v_c;
ierr = MatSetValuesStencil(jac, 1, &row, ncols, col, v, INSERT_VALUES); CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
MatNullSpace nullspace;
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);CHKERRQ(ierr);
ierr = MatSetNullSpace(jac,nullspace); CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&nullspace); CHKERRQ(ierr);
return 0;
}
#undef __FUNCT__
#define __FUNCT__ "setup"
PetscErrorCode setup(UserContext * ctx) {
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr);
ierr = KSPSetDMActive(ksp, PETSC_FALSE); CHKERRQ(ierr);
ierr = KSPSetDM(ksp, da_sc); CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp, ComputeMatrix, ctx); CHKERRQ(ierr);
ierr = KSPSetComputeRHS(ksp, ComputeRHS, ctx); CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
ierr = KSPSetUp(ksp); CHKERRQ(ierr);
return 0;
}
int main(int argc, char **argv) {
PetscInitialize(&argc, &argv, (char*)0, help);
UserContext ctx;
Vec phi;
PetscInt mx, my;
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
++rank;
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR,
128, 128, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da_sc); CHKERRQ(ierr);
ierr = DMDASetUniformCoordinates(da_sc, 0., 1., 0., 1., 0., 1.); CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da_sc, &ctx.rhs); CHKERRQ(ierr);
ierr = DMDAGetInfo(da_sc, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr);
ierr = VecSet(ctx.rhs, 1); CHKERRQ(ierr);
ierr = setup(&ctx); CHKERRQ(ierr);
PetscViewer viewer;
PetscViewerDrawOpen(PETSC_COMM_WORLD,NULL,NULL,0,0,300,300,&viewer);
for(int i = 0; i < N; ++i) {
ierr = VecSet(ctx.rhs, (float) i * mx); CHKERRQ(ierr);
ierr = VecScale(ctx.rhs, rank); CHKERRQ(ierr);
ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
ierr = KSPGetSolution(ksp, &phi); CHKERRQ(ierr);
ierr = VecView(ctx.rhs, viewer); CHKERRQ(ierr);
ierr = VecView(phi, viewer); CHKERRQ(ierr);
}
PetscViewerDestroy(&viewer);
VecDestroy(&phi);
PetscFinalize();
return 0;
}
