gistfile2.txt

sqrt_fib.cpp

#include "petscda.h"
#include "petscsnes.h"

typedef struct {
  DA da;
} UserCtx;

PetscErrorCode my_function(SNES snes,Vec x,Vec f,void *ctx);

#undef __FUNCT__  
#define __FUNCT__ "main"
int main(int argc, char *argv[])
{
  PetscErrorCode ierr;
  SNES snes;
  Mat J;
  Vec f, x;
  UserCtx *user;

  ierr = PetscInitialize(&argc, &argv, (char *)0, "Testing PETSc");CHKERRQ(ierr);
  ierr = PetscMalloc(sizeof(UserCtx),&user);CHKERRQ(ierr);
  ierr = DACreate1d(PETSC_COMM_WORLD, DA_XPERIODIC,-20,1,1,PETSC_NULL, &user->da);CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);

  ierr = DACreateGlobalVector(user->da, &x);CHKERRQ(ierr);
  ierr = VecDuplicate(x, &f);CHKERRQ(ierr);

  ierr = SNESSetFunction(snes, f, my_function, user);CHKERRQ(ierr);

  ierr = MatCreateSNESMF(snes, &J);CHKERRQ(ierr);
  ierr = SNESSetJacobian(snes, J, J, MatMFFDComputeJacobian, PETSC_NULL);CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  ierr = VecSet(x,1.0);CHKERRQ(ierr);
  ierr = SNESSolve(snes,0,x);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = MatDestroy(J);CHKERRQ(ierr);
  ierr = VecDestroy(f);CHKERRQ(ierr);
  ierr = VecDestroy(x);CHKERRQ(ierr);
  ierr = SNESDestroy(snes);CHKERRQ(ierr);
  ierr = PetscFree(user);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

#undef __FUNCT__  
#define __FUNCT__ "my_function"
PetscErrorCode my_function(SNES snes,Vec x,Vec f,void *ctx)
{
  PetscErrorCode ierr;
  UserCtx *user = (UserCtx*)ctx;
  DA da = user->da;
  Vec xloc;
  double *xa, *fa;
  int i0, i1, n, NX;

  PetscFunctionBegin;
  ierr = DAGetInfo(da,0, &NX,0,0, 0,0,0, 0,0,0,0);CHKERRQ(ierr);
  ierr = DAGetCorners(da, &i0,0,0, &n,0,0);CHKERRQ(ierr);
  i1 = i0 + n;

  ierr = DAGetLocalVector(da,&xloc);CHKERRQ(ierr);
  ierr = DAGlobalToLocalBegin(da,x,INSERT_VALUES,xloc);CHKERRQ(ierr);
  ierr = DAGlobalToLocalEnd(da,x,INSERT_VALUES,xloc);CHKERRQ(ierr);
  ierr = DAVecGetArray(da, xloc, &xa);CHKERRQ(ierr);
  ierr = DAVecGetArray(da, f, &fa);CHKERRQ(ierr);
//   double xb1 = 1.0*NX/(2*NX+1);
//   double xb2 = 1.0*NX/(2*NX+2);
//   fa[0] = NX*(xb2 - 6*xb1 + 3*xa[0] + 2*xa[1])/6 - xa[0]*xa[0];
//   fa[1] = NX*(xb1 - 6*xa[0] + 3*xa[1] + 2*xa[2])/6 - xa[1]*xa[1];
//   for (int i = 2; i < NX-1; i++) {
//     fa[i] = NX*(xa[i-2] - 6*xa[i-1] + 3*xa[i] + 2*xa[i+1])/6 - xa[i]*xa[i];
//   }
//   fa[NX-1] = NX*(xa[NX-3] - 6*xa[NX-2] + 3*xa[NX-1] + 2)/6 - xa[NX-1]*xa[NX-1];

// x = harmonic mean of neighbours

//   if (i0 == 0) {
//     fa[0] = xa[0]*xa[0] - 0.5*(xa[1]*xa[1]);
//     ii0 = 1;
//   } else if (i1 == NX) {
//     fa[NX-1] = xa[NX-1]*xa[NX-1] - 0.5*(xa[NX-2]*xa[NX-2] + 4.0);
//     ii1 = NX-1;
//   }
//   for (int i = ii0; i < ii1; i++) {
//     fa[i] = xa[i]*xa[i] - 0.5*(xa[i-1]*xa[i-1] + xa[i+1]*xa[i+1]);
//   }

// x = arithmetic mean of neighbours

  if (i0 == 0) {
    xa[-1] = 0;
  }
  if (i1 == NX) {
    xa[NX] = 1;
  }
  for (int i = i0; i < i1; i++) {
    fa[i] = xa[i] - 0.5*(xa[i-1] + xa[i+1]);
  }

  ierr = DAVecRestoreArray(da, f, &fa);CHKERRQ(ierr);
  ierr = DAVecRestoreArray(da, x, &xa);CHKERRQ(ierr);
  ierr = DARestoreLocalVector(da,&xloc);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}