On Thu, Jul 12, 2018 at 6:30 AM Ali Reza Khaz'ali
<[email protected] <mailto:[email protected]>> wrote:
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.
Why not just use DMDA? There are lots of examples (SNES ex5 for
instance) and it
handles the parallelism for you. This problem could be done in a few
hours I think.
Thanks,
Matt
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;
}
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/%7Emk51/>
