[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
code, please.
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 = VecGetArrayRead(x, &xx); CHKERRQ(ierr);
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 = VecRestoreArrayRead(x, &xx); CHKERRQ(ierr);
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
*/
ierr = VecGetArrayRead(x, &xx); CHKERRQ(ierr);
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
*/
ierr = VecRestoreArrayRead(x, &xx); CHKERRQ(ierr);
/*
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;
}
More information about the petsc-users
mailing list