Hello all,
I'm a relatively new PETSc user and am currently attempting to apply
matrix-free methods to a program I wrote that uses the KSP and DMDA
libraries.
However, I've come across an issue. When I use KSPSetDM, it requires that I
set operators with KSPSetComputeOperators rather than simply setting an
operator. However, the primary issue here is that when I use
KSPSetComputeOperators, I attempt to compute my operators with the
following function:
PetscErrorCode ComputeMatrix(KSP ksp, Mat A, Mat jac, void *user){
PetscErrorCode ierr;
AppCtx *ctx = (AppCtx*)user;
PetscFunctionBegin;
ierr = MatSetType(A, MATSHELL);CHKERRQ(ierr);
ierr =
MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &A);CHKERRQ(ierr);
ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MyMatMult);
CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = MatSetType(jac, MATSHELL);CHKERRQ(ierr);
ierr =
MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &jac);CHKERRQ(ierr);
ierr = MatShellSetOperation(jac,MATOP_MULT,(void(*)(void))MyMatMult);
CHKERRQ(ierr);
ierr = MatSetUp(jac);CHKERRQ(ierr);
PetscFunctionReturn(ierr);
}
When I try to run KSPSolve, I get the following error message:
joseph@Otto ~/petsc-3.9.3/Diffusion $ mpirun -n 1 ./diffu2matfree -x 3 -y 3
shell
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Must call MatXXXSetPreallocation() or MatSetUp() on
argument 1 "mat" before MatMult()
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.9.3, Jul, 02, 2018
[0]PETSC ERROR: ./diffu2matfree on a arch-linux2-c-debug named Otto by
joseph Thu Jul 19 17:53:00 2018
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
--with-fc=gfortran --download-mpich --download-fblaslapack
[0]PETSC ERROR: #1 MatMult() line 2306 in
/home/joseph/petsc-3.9.3/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 PCApplyBAorAB() line 682 in
/home/joseph/petsc-3.9.3/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #3 KSP_PCApplyBAorAB() line 304 in
/home/joseph/petsc-3.9.3/include/petsc/private/kspimpl.h
[0]PETSC ERROR: #4 KSPGMRESCycle() line 152 in
/home/joseph/petsc-3.9.3/src/ksp/ksp/impls/gmres/gmres.c
[0]PETSC ERROR: #5 KSPSolve_GMRES() line 234 in
/home/joseph/petsc-3.9.3/src/ksp/ksp/impls/gmres/gmres.c
[0]PETSC ERROR: #6 KSPSolve() line 669 in
/home/joseph/petsc-3.9.3/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #7 main() line 110 in
/home/joseph/petsc-3.9.3/Diffusion/diffu2matfree.cpp
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -x 3
[0]PETSC ERROR: -y 3
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to [email protected]
application called MPI_Abort(MPI_COMM_WORLD, 73) - process 0
Does anybody know what's going on here? I've also tried to make the
algorithm work without KSPSetDM and KSPSetComputeOperators, but my MatMult
step is dependent on DMDA, which is no longer used for the vectors in the
KSP if I don't use KSPSetDM. I've attached my code in full if there are any
questions as to the structure of my code.
Thanks,
Joseph Essman
static char help[] = "Solves 2d diffusion equation using distributed arrays \
as well as Krylov Subspace Methods. \n\n\
Parameters:\n\
-a: diffusion constant. Set to 30000 by default. \n\
-x: the x dimension of the space in which diffusion is taking place. Must be \
an integer value. Set to 2000 by default. \n\
-y: the y dimension of the space in which diffusion is taking place. Must be \
an integer value. Set to 100 by default. \n\
-delxy: the spatial interval between grid points. Set to 0.5 by default. \n\
-endtime: the total time over which the diffusion equation is solved. Set to \
0.2 by default.\n\
-delt: the time step size for solving the diffusion equation. Set to 0.001 \
by default.\n\
-interval: the interval between time steps written to files. Must be an integer\
value. Set to 3 by default.\n\n";
#include <petscdm.h>
#include <petscdmda.h>
#include <petscksp.h>
#include <string>
#include <sstream>
typedef struct {
PetscInt nx, ny, nt, size, loywid, loxwid, localsize, interval, low, high;
PetscReal f, delt;
Vec local, alsolocal, global;
} AppCtx;
extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void*);
extern PetscErrorCode ComputeInitCon(DM, Vec, void*);
int WriteGrid(AppCtx*, PetscInt, Vec);
PetscErrorCode MyMatMult(Mat,Vec,Vec);
int main(int argc, char **argv){
//Initialize function---------------------------------------------------------
//Declare variables
PetscErrorCode ierr;
KSP ksp;
PC pc;
AppCtx ctx;
PetscReal h = 0.5, delt = 0.001, endt = 0.2, alph = 30000;
PetscInt xlen = 2000, ylen = 100, i, intvl = 3;
DMBoundaryType bo = DM_BOUNDARY_NONE;
DMDAStencilType stype = DMDA_STENCIL_BOX;
Vec x, y;
Mat A;
DM da;
//Initialize PETSc
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
//Import variables from command line
ierr = PetscOptionsGetReal(NULL,NULL,"-a",&alph,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-delxy",&h,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-endtime",&endt,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-delt",&delt,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-x",&xlen,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-y",&ylen,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-interval",&intvl,NULL);CHKERRQ(ierr);
//Assign values to variables in application context
ctx.nx = (xlen/h) + 1; ctx.ny = (ylen/h) + 1; ctx.nt = (endt / delt) + 1;
ctx.size = ctx.nx * ctx.ny; ctx.f = (alph*delt) / (h*h); ctx.interval = intvl;
ctx.delt = delt;
//Prepare Krylov method solver (KSP) and distributed array (DMDA)-------------
//Create KSP and DMDA contexts
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = DMDACreate2d(PETSC_COMM_WORLD,bo,bo,stype,ctx.nx,ctx.ny,PETSC_DECIDE,\
PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);
//Configure DMDA
ierr = DMSetFromOptions(da);CHKERRQ(ierr);
ierr = DMSetUp(da);CHKERRQ(ierr);
ierr = DMDASetUniformCoordinates(da,0,xlen,0,ylen,0,0);CHKERRQ(ierr);
//Create solution vector
ierr = DMCreateGlobalVector(da, &x);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da, &ctx.global);CHKERRQ(ierr);
ierr = DMCreateLocalVector(da, &ctx.local);CHKERRQ(ierr);
ierr = DMCreateLocalVector(da, &ctx.alsolocal);CHKERRQ(ierr);
//Get local ownership range for local computations
ierr = VecGetOwnershipRange(x, &ctx.low, &ctx.high);CHKERRQ(ierr);
//Compute initial condition
ierr = ComputeInitCon(da, x, &ctx);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp, ComputeMatrix, &ctx);CHKERRQ(ierr);
//Set preconditioning (PC) method and configure KSP
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetUp(ksp);CHKERRQ(ierr);
ierr = KSPGetOperators(ksp, 0, &A);CHKERRQ(ierr);
MatType type;
MatGetType(A, &type);
PetscPrintf(PETSC_COMM_WORLD, "%s\n", type);
//Write initial condition to file
ierr = WriteGrid(&ctx, 0, x);CHKERRQ(ierr);
//Solve 2d diffusion equation-------------------------------------------------
//Iterate through all time steps
for (i=1; i<ctx.nt; i++){
//Solve at timestep
ierr = KSPSolve(ksp,x,x);CHKERRQ(ierr);
//Write solution to file
if(i%intvl == 0){
ierr = WriteGrid(&ctx, i, x);CHKERRQ(ierr);
}
}
//Wrap up function------------------------------------------------------------
//Clear memory
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
//Finalilze PETSc and return ierr
ierr = PetscFinalize();
return ierr;
}
PetscErrorCode ComputeInitCon(DM da, Vec b, void *ctx){
//Computes the vector that serves as the initial right hand side for the
//linear solver. This serves as our initial condition. In this case, the
//initial condition is a rectangle in the middle of the mesh, taking up about
//a quarter of the available space.
//Initialize function---------------------------------------------------------
//Declare variables
PetscErrorCode ierr;
PetscInt i,j,nx,ny,xm,ym,xs,ys,xlow,xhigh,ylow,yhigh;
PetscScalar **array;
//Begin function
PetscFunctionBegin;
//Prepare values and array for vector computation-----------------------------
//Create boundaries for rectangle
ierr = DMDAGetInfo(da, 0, &nx, &ny, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
xlow = (nx-1)/4; xhigh = nx - xlow-1; if(xhigh==nx){xhigh = nx - 1;}
ylow = (ny-1)/4; yhigh = ny - ylow-1; if(yhigh==ny){yhigh = ny - 1;}
//Get local boundaries
ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
//Compute initial condition vector--------------------------------------------
//Prepare array for value allocation
ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
//Iterate through local grid and populate array according to location
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
if (i>xlow && i<xhigh && j>ylow &&j<yhigh){array[j][i] = 1;}
else {array[j][i] = 0;}
}
}
//Restore array to vector and assemble completed right hand side
ierr = DMDAVecRestoreArray(da, b, &array);CHKERRQ(ierr);
ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
//Return error code
PetscFunctionReturn(ierr);
}
PetscErrorCode ComputeMatrix(KSP ksp, Mat A, Mat jac, void *user){
PetscErrorCode ierr;
AppCtx *ctx = (AppCtx*)user;
PetscFunctionBegin;
ierr = MatSetType(A, MATSHELL);CHKERRQ(ierr);
ierr = MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &A);CHKERRQ(ierr);
ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MyMatMult);
CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = MatSetType(jac, MATSHELL);CHKERRQ(ierr);
ierr = MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &jac);CHKERRQ(ierr);
ierr = MatShellSetOperation(jac,MATOP_MULT,(void(*)(void))MyMatMult);
CHKERRQ(ierr);
ierr = MatSetUp(jac);CHKERRQ(ierr);
PetscFunctionReturn(ierr);
}
int WriteGrid(AppCtx *ctx, PetscInt n, Vec u){
//Writes distributed array grid to a binary file readable by MATLAB
//Initialize function---------------------------------------------------------
//Declare variables
PetscErrorCode ierr;
PetscViewer viewer;
char filename[100], num[15];
//Begin function
PetscFunctionBegin;
//Write main grid-------------------------------------------------------------
//Prepare file name
strcpy(filename, "VTKFolder/grid");
sprintf(num, "%d.vtr", n);
strcat(filename, num);
//Open file and write vector
ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTR);CHKERRQ(ierr);
ierr = VecView(u,viewer);CHKERRQ(ierr);
//Free memory and return error code-------------------------------------------
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
PetscFunctionReturn(ierr);
}
PetscErrorCode MyMatMult(Mat A, Vec X, Vec Y){
PetscErrorCode ierr;
void *ptr;
AppCtx *user;
PetscScalar **yarray, **xarray;
DM da;
PetscInt i,j,xm,ym,xs,ys,gxm,gym,gxs,gys;
PetscFunctionBegin;
ierr = MatShellGetContext(A,&ptr);CHKERRQ(ierr);
ierr = VecGetDM(X,&da);CHKERRQ(ierr);
user = (AppCtx*)ptr;
PetscScalar f = user->f;
ierr = DMGlobalToLocalBegin(da, X, INSERT_VALUES, user->local);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(da, X, INSERT_VALUES, user->local);CHKERRQ(ierr);
ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gxm,&gym,0);CHKERRQ(ierr);
ierr = DMDAVecGetArray(da, Y, &yarray);CHKERRQ(ierr);
ierr = DMDAVecGetArrayRead(da, user->local, &xarray);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
if (i == gxs || j == gys || i == gxs+gxm-1 || j == gys+gym-1){
yarray[j][i] = xarray[j][i];
}
else{
yarray[j][i] = -f*xarray[j-1][i]-f*xarray[j+1][i]-f*xarray[j][i-1]\
-f*xarray[j][i+1]+(1+4*f)*xarray[j][i];
}
}
}
ierr = DMDAVecRestoreArray(da, Y, &yarray);CHKERRQ(ierr);
ierr = DMDAVecRestoreArrayRead(da, user->local, &xarray);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
