# [petsc-users] Problem in SNES

Matthew Knepley knepley at gmail.com
Thu Jul 12 10:30:15 CDT 2018

```On Thu, Jul 12, 2018 at 11:14 AM Ali Reza Khaz'ali <arkhazali at cc.iut.ac.ir>
wrote:

>
> I am experimenting with SNES for now, however, my purpose is to solve a
> multi-physics large system with unstructured grids (Enhanced Oil Recovery
> simulation in fractured petroleum reservoirs to be specific, phase
> equilibrium and three phase flow equations must be solved simultaneously),
> where to the extend of my knowledge, no existing PETSc vector structure can
> be used. Therefore, I started with solving simple ODEs with SNES without
> using those structures.
>
The intention is that you use DMPlex for this.

> An interesting fact about my code samples is, with increasing the number
> of MPI workers, the number of SNES iterations (and hence, the precision of
> the solution) decrease. I'm sure there is a flaw in my code (I do not have
> much MPI experience), but I can't find it.
>
We really do not have the resources to debug raw MPI code. That is why we
have the examples.

Thanks,

Matt

> Best wishes,
>
> Ali
>
> On 7/12/2018 3:58 PM, Matthew Knepley wrote:
>
> 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
>>
>
> 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
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/%7Emk51/>
>
>
>

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their