[petsc-users] Parallelizing Nonlinear Algebraic systems without DM

[Nidish]

Hello,

My question is related to solving a generic nonlinear algebraic system 
in parallel. I'm using the word generic to point to the fact that the 
system may or may not be from a finite element/mesh type object - it's 
just a residual function and a Jacobian function that I can provide.

I followed the serial example in $PETSC_DIR/src/snes/tutorials/ex1.c 
pretty much verbatim (albeit with a slightly different objective 
function whose size can be controlled), and have been trying to extend 
it to a parallel application (the given example is strictly serial).

When I run my code in serial everything works well and I get the 
solution I expect. However, even if I use just two cores, the code 
starts returning garbage. I had a feeling it has something to do with 
index ordering, but I want to get this code working without using a DM 
object. I'm attaching my code with this email.

TL;DR Version: What modifications must be made in 
"$PETSC_DIR/src/snes/tutorials/ex1.c" to make it work in Parallel? Is it 
possible to make minimal modifications to make it run in parallel too?

Thank You,
Nidish
static char help[] = "Nonlinear Solution Trial\n";

#include <petscsnes.h>
#include <petscdmda.h>

PetscErrorCode NLFUN(SNES, Vec, Vec, void*);
PetscErrorCode NLJAC(SNES, Vec, Mat, Mat, void*);

int main(int nargs, char *sargs[])
{
  PetscInt ierr;
  PetscMPIInt rank, size;
  ierr = PetscInitialize(&nargs, &sargs, (char*)0, help); if (ierr) return ierr;
  MPI_Comm_size(PETSC_COMM_WORLD, &size);
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

  PetscInt N=5;
  PetscScalar Fb[] = {30, 2.5};
  SNES snes;
  KSP ksp;
  PC pc;
  Vec x, r;
  Mat J;
  PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
  PetscOptionsGetScalar(NULL, NULL, "-F", &Fb[0], NULL);
  PetscOptionsGetScalar(NULL, NULL, "-b", &Fb[1], NULL);

  PetscPrintf(PETSC_COMM_WORLD, "%d\n", N);
  
  SNESCreate(PETSC_COMM_WORLD, &snes);

  MatCreate(PETSC_COMM_WORLD, &J);
  MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, N, N);
  MatSetFromOptions(J);
  MatSetUp(J);

  MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
  
  MatCreateVecs(J, &x, &r);
  VecSetFromOptions(x);
  VecSetFromOptions(r);

  /* PetscInt m, n; */
  /* MatGetLocalSize(J, &m, &n); */
  /* MatMPIAIJSetPreallocation(J, 3, NULL, N-3, NULL); */
  
  SNESSetFunction(snes, r, NLFUN, (void*)Fb);
  SNESSetJacobian(snes, J, J, NLJAC, (void*)Fb);

  SNESGetKSP(snes, &ksp);
  KSPGetPC(ksp, &pc);
  PCSetType(pc, PCNONE);
  PCSetFromOptions(pc);
  KSPSetFromOptions(ksp);

  SNESSetFromOptions(snes);

  /* Initial Guess ***********************************************************/
  /* VecSetValue(r, N-1, Fb, INSERT_MODE); */
  VecSet(r, 0.0);
  
  VecAssemblyBegin(r);
  
  VecSet(x, sqrt((PetscScalar)N));

  VecAssemblyBegin(x);
  VecAssemblyEnd(r);
  VecAssemblyEnd(x);

  ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);

  VecAssemblyBegin(x);
  VecAssemblyBegin(r);
  VecAssemblyEnd(x);
  VecAssemblyEnd(r);

  SNESDestroy(&snes);
  VecDestroy(&x);
  MatDestroy(&J);
  PetscFinalize();
  return 0;
}

PetscErrorCode NLFUN(SNES snes, Vec x, Vec f, void* ctx)
{
  PetscErrorCode ierr;
  /* const PetscScalar *xx; */
  PetscScalar* xx;
  PetscScalar* Fb = (PetscScalar*) ctx;
  /* PetscScalar xnrm; */
  
  PetscInt is, ie, N;
  VecGetSize(x, &N);
  VecGetOwnershipRange(x, &is, &ie);
  /* VecGetArrayRead(x, &xx); */
  VecGetArray(x, &xx);  
  for (int i=is; i<ie; ++i)
    if (i==0)
      VecSetValue(f, i, 2*xx[i -is]-xx[i+1 -is], INSERT_VALUES);
    else if (i<N-1)
      VecSetValue(f, i, -xx[i-1 -is]+2*xx[i -is]-xx[i+1 -is], INSERT_VALUES);
    else if (i==N-1) {
      /* VecNorm(x, NORM_2, &xnrm); */
      /* VecSetValue(f, i, 0.5*(xnrm*xnrm-N*N), INSERT_VALUES); */

      VecSetValue(f, i, -xx[i-1 -is]+xx[i -is]+Fb[1]*pow(xx[i -is], 3)-Fb[0], INSERT_VALUES);
    }
  VecAssemblyBegin(f);
  /* VecRestoreArrayRead(x, &xx); */
  VecRestoreArray(x, &xx);  
  ierr = VecAssemblyEnd(f);

  return ierr;
}

PetscErrorCode NLJAC(SNES snes, Vec x, Mat jac, Mat B, void* ctx)
{
  PetscErrorCode ierr;
  PetscScalar xxend;
  PetscScalar* Fb = (PetscScalar*) ctx;  
  /* PetscScalar xnrm; */

  PetscInt is, ie, N;
  VecGetSize(x, &N);
  VecGetOwnershipRange(x, &is, &ie);

  PetscInt Nend = N-1;
  for (int i=is; i<ie; ++i)
    if (i==0) {
      MatSetValue(B, i, i, 2.0, INSERT_VALUES);
      MatSetValue(B, i, i+1, -1.0, INSERT_VALUES);      
    }
    else if (i<N-1) {
      MatSetValue(B, i, i-1, -1.0, INSERT_VALUES);
      MatSetValue(B, i, i, 2.0, INSERT_VALUES);
      MatSetValue(B, i, i+1, -1.0, INSERT_VALUES);
    }
    else if (i==N-1) {
      /* VecNorm(x, NORM_2, &xnrm); */
      /* MatSetValuesRow(jac, i, xx); */
      
      VecGetValues(x, 1, &Nend, &xxend);
      
      MatSetValue(B, i, i-1, -1.0, INSERT_VALUES);
      MatSetValue(B, i, i, 1.0+3*Fb[1]*pow(xxend, 2), INSERT_VALUES);

      /* PetscPrintf(PETSC_COMM_WORLD, "%e\n", xxend); */
    }
  MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
  ierr = MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

  if (jac != B) {
  MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);    
  }
  
  return ierr;
}