# [petsc-users] Problem in SNES

Ali Reza Khaz'ali arkhazali at cc.iut.ac.ir
Thu Jul 12 05:29:57 CDT 2018

```Hi,

I have to apologize if I'm asking a noob question. I am trying to solve
y"-(1-x/5)yy'=x with boundary condition of y(1)=2 and y(3)=-1 with
finite difference method. The discretized equations is to be solved by
SNES. I wrote the follwoing code, however, it only works if number of
MPI processes is set to one. for more workers, it does not seem to give
valid results. I'd be grateful if anybody can guide me to debug the

Regards,

Ali

static char help[] = "Newton's method for a non-linear ODE system.\n\n";
double h;
int nParts;

#include <petscsnes.h>

/*
User-defined routines
*/
extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void*);
extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);

int main(int argc, char **argv)
{
SNES           snes;         /* nonlinear solver context */
KSP            ksp;          /* linear solver context */
PC             pc;           /* preconditioner context */
Vec            x, r;          /* solution, residual vectors */
Mat            J;            /* Jacobian matrix */
PetscErrorCode ierr;
PetscMPIInt    size;
PetscScalar    pfive = 0, *xx;
int rank;
int i, j;
Vec f;

nParts = 30;
h = (3.0 - 1.0) / nParts;

ierr = PetscInitialize(&argc, &argv, (char*)0, help); if (ierr)
return ierr;
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
ierr = SNESCreate(PETSC_COMM_WORLD, &snes); CHKERRQ(ierr);
ierr = VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, nParts-1, &x);
CHKERRQ(ierr);
ierr = VecDuplicate(x, &r); CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,
nParts - 1, nParts - 1, 3, NULL, 3, NULL, &J); CHKERRQ(ierr);

ierr = SNESSetFunction(snes, r, FormFunction, NULL); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, J, J, FormJacobian, NULL); CHKERRQ(ierr);

ierr = SNESGetKSP(snes, &ksp); CHKERRQ(ierr);
ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp, 1.e-8, 1e-16, PETSC_DEFAULT, 200);
CHKERRQ(ierr);
ierr = SNESSetTolerances(snes, 1.e-8, 1e-16, PETSC_DEFAULT, 200,
200); CHKERRQ(ierr);
ierr = SNESSetType(snes, SNESNEWTONLS); CHKERRQ(ierr);
ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr);

ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
ierr = VecSet(x, pfive); CHKERRQ(ierr);
ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);
ierr = SNESGetFunction(snes, &f, 0, 0); CHKERRQ(ierr);

ierr = VecDestroy(&x); CHKERRQ(ierr); ierr = VecDestroy(&r);
CHKERRQ(ierr);
ierr = MatDestroy(&J); CHKERRQ(ierr); ierr = SNESDestroy(&snes);
CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}

PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
{
PetscErrorCode    ierr;
const PetscScalar *xx;
PetscScalar       *ff;
PetscInt vstart, vend;
int i, j, rank, size, xlength;
double xp, xxstart, xxend, x1, x2;
MPI_Status MyStat;

xxstart = 0;
xxend = 0;

MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);

ierr = VecGetOwnershipRange(x, &vstart, &vend); CHKERRQ(ierr);
xlength = vend - vstart;
ierr = VecGetArray(f, &ff); CHKERRQ(ierr);
ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD);
if (size > 1) {
if (rank == 0) {
ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, 1, 0,
PETSC_COMM_WORLD);
ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, 1, 1,
PETSC_COMM_WORLD, &MyStat);
}
else if (rank == (size - 1)) {
ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1, rank -
1, PETSC_COMM_WORLD, &MyStat);
ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
PETSC_COMM_WORLD);
}
else {
ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, rank + 1,
rank, PETSC_COMM_WORLD);
ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, rank + 1, rank + 1,
PETSC_COMM_WORLD, &MyStat);
ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
PETSC_COMM_WORLD);
ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1, rank -
1, PETSC_COMM_WORLD, &MyStat);
}
}

for (i = vstart; i < vend; i++) {
xp = 1 + i * h;
j = i - vstart;

if (i == vstart) x1 = xxstart;
else x1 = xx[j - 1];
if (i == (vend - 1)) x2 = xxend;
else x2 = xx[j + 1];

if (i==0) ff[j] = (x2 - 2 * xx[j] + 2) / (h*h) - (1 - xp /
5)*xx[j] * (x2 - 2) / (2 * h) - xp;
else if (i==(nParts-2)) ff[j]=(-1 - 2 * xx[j] + x1) / (h*h) -
(1 - xp / 5)*xx[j] * (-1 - x1) / (2 * h) - xp;
else ff[j] = (x2 - 2 * xx[j] + x1) / (h*h) - (1 - xp / 5)*xx[j]
* (x2 - x1) / (2 * h) - xp;
}

/* Restore vectors */
ierr = VecRestoreArray(f, &ff); CHKERRQ(ierr);
return 0;
}

PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
{
const PetscScalar *xx;
PetscErrorCode    ierr;

PetscInt vstart, vend;
double rrow[3];
int col[3];
int i, j, rank, size, xlength;
double xp, xxstart, xxend, x1, x2;
MPI_Status MyStat;

ierr = MatZeroEntries(jac); CHKERRQ(ierr);

ierr = VecGetOwnershipRange(x, &vstart, &vend); CHKERRQ(ierr);
xlength = vend - vstart;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);

xxstart = 0;
xxend = 0;

/*
Get pointer to vector data.0+1
*/
if (size > 1) {
if (rank == 0) {
ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, 1, 0,
PETSC_COMM_WORLD);
ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, 1, 1,
PETSC_COMM_WORLD, &MyStat);
xxstart = 1e12;
}
else if (rank == (size - 1)) {
ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1, rank -
1, PETSC_COMM_WORLD, &MyStat);
ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
PETSC_COMM_WORLD);
xxend = 1e12;
}
else {
ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, rank + 1,
rank, PETSC_COMM_WORLD);
ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, rank + 1, rank + 1,
PETSC_COMM_WORLD, &MyStat);
ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
PETSC_COMM_WORLD);
ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1, rank -
1, PETSC_COMM_WORLD, &MyStat);
}
}

for (i = vstart; i < vend; i++) {
xp = 1 + i * h;
j = i - vstart;

if (i == vstart) x1 = xxstart;
else x1 = xx[j - 1];
if (i == (vend - 1)) x2 = xxend;
else x2 = xx[j + 1];

if (i == 0) {
rrow[0] = -2 / (h*h) - (1 - xp / 5)*(x2 - 2) / (2 * h);
rrow[1] = 1 / (h*h) - (1 + xp / 5)*xx[j] / (2 * h);
col[0] = i;
col[1] = i + 1;
ierr = MatSetValues(jac, 1, &i, 2, col, rrow,
INSERT_VALUES); CHKERRQ(ierr);
}
else if (i == (nParts-2)) {
rrow[0] = 1 / (h*h) + (1 + xp / 5)*xx[j] / (2 * h);
rrow[1] = -2 / (h*h) - (1 - xp / 5)*(-1 - x1) / (2 * h);
col[0] = i - 1;
col[1] = i;
ierr = MatSetValues(jac, 1, &i, 2, col, rrow,
INSERT_VALUES); CHKERRQ(ierr);
}
else {
rrow[0] = 1 / (h*h) + (1 + xp / 5)*xx[j] / (2 * h);
rrow[1] = -2 / (h*h) - (1 - xp / 5)*(x2 - x1) / (2 * h);
rrow[2] = 1 / (h*h) - (1 + xp / 5)*xx[j] / (2 * h);
col[0] = i-1;
col[1] = i;
col[2] = i + 1;
ierr = MatSetValues(jac, 1, &i, 3, col, rrow,
INSERT_VALUES); CHKERRQ(ierr);
}
}

/*
Restore vector
*/

/*
Assemble matrix
*/
B = jac;
ierr = MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

return 0;
}

```