[petsc-users] Problem in SNES
Matthew Knepley
knepley at gmail.com
Thu Jul 12 06:28:11 CDT 2018
On Thu, Jul 12, 2018 at 6:30 AM Ali Reza Khaz'ali <arkhazali at cc.iut.ac.ir>
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/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180712/c3a5fbb8/attachment.html>
More information about the petsc-users
mailing list