C++ SNESSolve函数代码示例

unsigned int PetscDiffSolver::solve()
{
  this->init();

  START_LOG("solve()", "PetscDiffSolver");

  PetscVector<Number> &x =
    *(libmesh_cast_ptr<PetscVector<Number>*>(_system.solution.get()));
  PetscMatrix<Number> &jac =
    *(libmesh_cast_ptr<PetscMatrix<Number>*>(_system.matrix));
  PetscVector<Number> &r =
    *(libmesh_cast_ptr<PetscVector<Number>*>(_system.rhs));

#ifdef LIBMESH_ENABLE_CONSTRAINTS
  _system.get_dof_map().enforce_constraints_exactly(_system);
#endif

  int ierr = 0;

  ierr = SNESSetFunction (_snes, r.vec(),
                          __libmesh_petsc_diff_solver_residual, this);
    LIBMESH_CHKERRABORT(ierr);

  ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(),
                          __libmesh_petsc_diff_solver_jacobian, this);
    LIBMESH_CHKERRABORT(ierr);

# if PETSC_VERSION_LESS_THAN(2,2,0)

  ierr = SNESSolve (_snes, x.vec(), &_outer_iterations);
         LIBMESH_CHKERRABORT(ierr);

// 2.2.x style
#elif PETSC_VERSION_LESS_THAN(2,3,0)

  ierr = SNESSolve (_snes, x.vec());
         LIBMESH_CHKERRABORT(ierr);

// 2.3.x & newer style
#else

  ierr = SNESSolve (_snes, PETSC_NULL, x.vec());
         LIBMESH_CHKERRABORT(ierr);

#endif

  STOP_LOG("solve()", "PetscDiffSolver");

  SNESConvergedReason reason;
  SNESGetConvergedReason(_snes, &reason);

  this->clear();

  return convert_solve_result(reason);
}

unsigned int PetscDiffSolver::solve()
{
  this->init();

  START_LOG("solve()", "PetscDiffSolver");

  PetscVector<Number> &x =
    *(libmesh_cast_ptr<PetscVector<Number>*>(_system.solution.get()));
  PetscMatrix<Number> &jac =
    *(libmesh_cast_ptr<PetscMatrix<Number>*>(_system.matrix));
  PetscVector<Number> &r =
    *(libmesh_cast_ptr<PetscVector<Number>*>(_system.rhs));

#ifdef LIBMESH_ENABLE_CONSTRAINTS
  _system.get_dof_map().enforce_constraints_exactly(_system);
#endif

  int ierr = 0;

  ierr = SNESSetFunction (_snes, r.vec(),
                          __libmesh_petsc_diff_solver_residual, this);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

  ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(),
                          __libmesh_petsc_diff_solver_jacobian, this);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

# if PETSC_VERSION_LESS_THAN(2,2,0)

  ierr = SNESSolve (_snes, x.vec(), &_outer_iterations);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.2.x style
#elif PETSC_VERSION_LESS_THAN(2,3,0)

  ierr = SNESSolve (_snes, x.vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.3.x & newer style
#else

  ierr = SNESSolve (_snes, PETSC_NULL, x.vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

#endif

  STOP_LOG("solve()", "PetscDiffSolver");

  this->clear();

  // FIXME - We'll worry about getting the solve result right later...

  return DiffSolver::CONVERGED_RELATIVE_RESIDUAL;
}

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    SNES      snes;          // nonlinear solver context
    Vec       x, r;          // solution, residual vectors
    PetscReal x0 = 2.0;

    PetscInitialize(&argc,&argv,(char*)0,help);

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","options to atan","");CHKERRQ(ierr);
    ierr = PetscOptionsReal("-x0","initial value","atan.c",x0,&x0,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd();CHKERRQ(ierr);

    ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,1); CHKERRQ(ierr);
    ierr = VecSetFromOptions(x); CHKERRQ(ierr);
    ierr = VecSet(x,x0); CHKERRQ(ierr);
    ierr = VecDuplicate(x,&r); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,r,FormFunction,NULL); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,x); CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

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

int main(int argc, char **argv)
{
  DM             dm;   /* Problem specification */
  SNES           snes; /* Nonlinear solver */
  Vec            u;    /* Solutions */
  AppCtx         user; /* User-defined work context */
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
  /* Primal system */
  ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
  ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
  ierr = SetupFE(dm, user.dim, user.simplex, "displacement", SetupPrimalProblem, &user);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
  ierr = VecSet(u, 0.0);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) u, "displacement");CHKERRQ(ierr);
  ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
  ierr = DMSNESCheckFromOptions(snes, u, NULL, NULL);CHKERRQ(ierr);
  ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
  ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr);
  ierr = VecViewFromOptions(u, NULL, "-displacement_view");CHKERRQ(ierr);
  /* Cleanup */
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

static PetscErrorCode TSStage_EIMEX(TS ts,PetscInt istage)
{
  TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
  PetscReal       h;
  Vec             Y=ext->Y, Z=ext->Z;
  SNES            snes;
  TSAdapt         adapt;
  PetscInt        i,its,lits;
  PetscBool       accept;
  PetscErrorCode  ierr;

  PetscFunctionBegin;
  ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
  h = ts->time_step/ext->N[istage];/* step size for the istage-th stage */
  ext->shift = 1./h;
  ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */
  ierr = VecCopy(ext->VecSolPrev,Y);CHKERRQ(ierr); /* Take the previous solution as intial step */

  for(i=0; i<ext->N[istage]; i++){
    ext->ctime = ts->ptime + h*i;
    ierr = VecCopy(Y,Z);CHKERRQ(ierr);/* Save the solution of the previous substep */
    ierr = SNESSolve(snes,NULL,Y);CHKERRQ(ierr);
    ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
    ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr);
    ts->snes_its += its; ts->ksp_its += lits;
    ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);
    ierr = TSAdaptCheckStage(adapt,ts,ext->ctime,Y,&accept);CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

// -------------------------------------------------------------
// PetscNonlinearSolverImplementation::p_solve
// -------------------------------------------------------------
void
PetscNonlinearSolverImplementation::p_solve(void)
{
  PetscErrorCode ierr(0);
  p_petsc_X = PETScVector(*p_X);
  int me(this->processor_rank());

  try {
    ierr = SNESSolve(p_snes, NULL, *p_petsc_X); CHKERRXX(ierr);
    SNESConvergedReason reason;
    PetscInt iter;
    ierr = SNESGetConvergedReason(p_snes, &reason); CHKERRXX(ierr);
    ierr = SNESGetIterationNumber(p_snes, &iter); CHKERRXX(ierr);

    std::string msg;
    if (reason < 0) {
      msg = 
        boost::str(boost::format("%d: PETSc SNES diverged after %d iterations, reason: %d") % 
                   me % iter % reason);
      throw Exception(msg);
    } else if (me == 0) {
      msg = 
        boost::str(boost::format("%d: PETSc SNES converged after %d iterations, reason: %d") % 
                   me % iter % reason);
      std::cerr << msg << std::endl;
    }
  } catch (const PETSC_EXCEPTION_TYPE& e) {
    throw PETScException(ierr, e);
  } catch (const Exception& e) {
    throw e;
  }
}

/*
Defines the action of the downsmoother
 */
PetscErrorCode SNESFASDownSmooth_Private(SNES snes, Vec B, Vec X, Vec F, PetscReal *fnorm)
{
  PetscErrorCode      ierr = 0;
  SNESConvergedReason reason;
  Vec                 FPC;
  SNES                smoothd;
  SNES_FAS            *fas = (SNES_FAS*) snes->data;

  PetscFunctionBegin;
  ierr = SNESFASCycleGetSmootherDown(snes, &smoothd);CHKERRQ(ierr);
  ierr = SNESSetInitialFunction(smoothd, F);CHKERRQ(ierr);
  if (fas->eventsmoothsolve) {ierr = PetscLogEventBegin(fas->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
  ierr = SNESSolve(smoothd, B, X);CHKERRQ(ierr);
  if (fas->eventsmoothsolve) {ierr = PetscLogEventEnd(fas->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
  /* check convergence reason for the smoother */
  ierr = SNESGetConvergedReason(smoothd,&reason);CHKERRQ(ierr);
  if (reason < 0 && !(reason == SNES_DIVERGED_MAX_IT || reason == SNES_DIVERGED_LOCAL_MIN || reason == SNES_DIVERGED_LINE_SEARCH)) {
    snes->reason = SNES_DIVERGED_INNER;
    PetscFunctionReturn(0);
  }
  ierr = SNESGetFunction(smoothd, &FPC, NULL, NULL);CHKERRQ(ierr);
  ierr = VecCopy(FPC, F);CHKERRQ(ierr);
  if (fnorm) {ierr = VecNorm(F,NORM_2,fnorm);CHKERRQ(ierr);}
  PetscFunctionReturn(0);
}

int main(int argc,char **argv)
{
  SNES           snes;
  PetscErrorCode ierr;
  PetscInt       its,lits;
  PetscReal      litspit;
  DM             da;

  PetscInitialize(&argc,&argv,PETSC_NULL,help);

  /*
      Set the DMDA (grid structure) for the grids.
  */
  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);
  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,PETSC_NULL);CHKERRQ(ierr);
  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);

  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  ierr = SNESSolve(snes,0,0);CHKERRQ(ierr);
  ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
  ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr);
  litspit = ((PetscReal)lits)/((PetscReal)its);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",litspit);CHKERRQ(ierr);

  ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  ierr = PetscFinalize();

  return 0;
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  UserCtx        user;
  DM             red,da;
  SNES           snes;
  DM             packer;
  PetscBool      use_monitor = PETSC_FALSE;

  ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
  ierr = PetscOptionsSetFromOptions(NULL);CHKERRQ(ierr);

  /* Hardwire several options; can be changed at command line */
  ierr = PetscOptionsInsertString(NULL,common_options);CHKERRQ(ierr);
  ierr = PetscOptionsInsertString(NULL,matrix_free_options);CHKERRQ(ierr);
  ierr = PetscOptionsInsert(NULL,&argc,&argv,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(NULL,NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);CHKERRQ(ierr);

  /* Create a global vector that includes a single redundant array and two da arrays */
  ierr = DMCompositeCreate(PETSC_COMM_WORLD,&packer);CHKERRQ(ierr);
  ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);CHKERRQ(ierr);
  ierr = DMSetOptionsPrefix(red,"red_");CHKERRQ(ierr);
  ierr = DMCompositeAddDM(packer,red);CHKERRQ(ierr);
  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-5,2,1,NULL,&da);CHKERRQ(ierr);
  ierr = DMSetOptionsPrefix(red,"da_");CHKERRQ(ierr);
  ierr = DMCompositeAddDM(packer,(DM)da);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(packer,&user);CHKERRQ(ierr);

  packer->ops->creatematrix = DMCreateMatrix_MF;

  /* create nonlinear multi-level solver */
  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes,packer);CHKERRQ(ierr);
  ierr = SNESSetFunction(snes,NULL,ComputeFunction,NULL);CHKERRQ(ierr);
  ierr = SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL);CHKERRQ(ierr);

  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  if (use_monitor) {
    /* create graphics windows */
    ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);CHKERRQ(ierr);
    ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);CHKERRQ(ierr);
    ierr = SNESMonitorSet(snes,Monitor,0,0);CHKERRQ(ierr);
  }

  ierr = SNESSolve(snes,NULL,NULL);CHKERRQ(ierr);
  ierr = SNESDestroy(&snes);CHKERRQ(ierr);

  ierr = DMDestroy(&red);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = DMDestroy(&packer);CHKERRQ(ierr);
  if (use_monitor) {
    ierr = PetscViewerDestroy(&user.u_lambda_viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&user.fu_lambda_viewer);CHKERRQ(ierr);
  }
  ierr = PetscFinalize();
  return ierr;
}

/*

Performs the FAS coarse correction as:

fine problem:   F(x) = b
coarse problem: F^c(x^c) = b^c

b^c = F^c(Rx) - R(F(x) - b)

 */
PetscErrorCode SNESFASCoarseCorrection(SNES snes, Vec X, Vec F, Vec X_new)
{
  PetscErrorCode      ierr;
  Vec                 X_c, Xo_c, F_c, B_c;
  SNESConvergedReason reason;
  SNES                next;
  Mat                 restrct, interpolate;
  SNES_FAS            *fasc;

  PetscFunctionBegin;
  ierr = SNESFASCycleGetCorrection(snes, &next);CHKERRQ(ierr);
  if (next) {
    fasc = (SNES_FAS*)next->data;

    ierr = SNESFASCycleGetRestriction(snes, &restrct);CHKERRQ(ierr);
    ierr = SNESFASCycleGetInterpolation(snes, &interpolate);CHKERRQ(ierr);

    X_c  = next->vec_sol;
    Xo_c = next->work[0];
    F_c  = next->vec_func;
    B_c  = next->vec_rhs;

    if (fasc->eventinterprestrict) {ierr = PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
    ierr = SNESFASRestrict(snes,X,Xo_c);CHKERRQ(ierr);
    /* restrict the defect: R(F(x) - b) */
    ierr = MatRestrict(restrct, F, B_c);CHKERRQ(ierr);
    if (fasc->eventinterprestrict) {ierr = PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}

    if (fasc->eventresidual) {ierr = PetscLogEventBegin(fasc->eventresidual,0,0,0,0);CHKERRQ(ierr);}
    /* F_c = F^c(Rx) - R(F(x) - b) since the second term was sitting in next->vec_rhs */
    ierr = SNESComputeFunction(next, Xo_c, F_c);CHKERRQ(ierr);
    if (fasc->eventresidual) {ierr = PetscLogEventEnd(fasc->eventresidual,0,0,0,0);CHKERRQ(ierr);}

    /* solve the coarse problem corresponding to F^c(x^c) = b^c = F^c(Rx) - R(F(x) - b) */
    ierr = VecCopy(B_c, X_c);CHKERRQ(ierr);
    ierr = VecCopy(F_c, B_c);CHKERRQ(ierr);
    ierr = VecCopy(X_c, F_c);CHKERRQ(ierr);
    /* set initial guess of the coarse problem to the projected fine solution */
    ierr = VecCopy(Xo_c, X_c);CHKERRQ(ierr);

    /* recurse to the next level */
    ierr = SNESSetInitialFunction(next, F_c);CHKERRQ(ierr);
    ierr = SNESSolve(next, B_c, X_c);CHKERRQ(ierr);
    ierr = SNESGetConvergedReason(next,&reason);CHKERRQ(ierr);
    if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
      snes->reason = SNES_DIVERGED_INNER;
      PetscFunctionReturn(0);
    }
    /* correct as x <- x + I(x^c - Rx)*/
    ierr = VecAXPY(X_c, -1.0, Xo_c);CHKERRQ(ierr);

    if (fasc->eventinterprestrict) {ierr = PetscLogEventBegin(fasc->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
    ierr = MatInterpolateAdd(interpolate, X_c, X, X_new);CHKERRQ(ierr);
    if (fasc->eventinterprestrict) {ierr = PetscLogEventEnd(fasc->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
  }
  PetscFunctionReturn(0);
}

/**
   In theory we have SolveTimeDependent and SolveSteadyState, the latter doesn't work very well.

   @param ctx the application context
*/
PetscErrorCode SolveTimeDependent(void* ctx)
{
  User                 user = (User)ctx;
  PetscMPIInt          rank, size;
  PetscErrorCode       ierr;

  Algebra              algebra = user->algebra;

  PetscLogDouble v1, v2;
  SNESConvergedReason snesreason;

  PetscFunctionBegin;

  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);

  /* set current time and step */
  user->current_time = user->initial_time;
  user->current_step = 1;
  ierr = SNESMonitorSet(user->snes, MonitorFunction, (void*) user,
			PETSC_NULL);CHKERRQ(ierr);

  /* start the time step iteration */
  while (user->current_time < (user->final_time - 0.05 * user->dt)) {
    /*Do a simple adaptive time step size*/

    ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);
    ierr = VecCopy(algebra->fn, algebra->oldfn);CHKERRQ(ierr);

    /* b) update the current time */
    user->current_time += user->dt;

    ierr = PetscPrintf(PETSC_COMM_WORLD,"\nCurrent time = %f\n",
		         user->current_time); CHKERRQ(ierr);

    /* c) solve the time dependent problem */
    ierr = PetscTime(&v1);CHKERRQ(ierr);
    ierr = SNESSolve(user->snes, PETSC_NULL, algebra->solution);CHKERRQ(ierr);
    ierr = PetscTime(&v2);CHKERRQ(ierr);

    ierr = SNESGetConvergedReason(user->snes,&snesreason);CHKERRQ(ierr);
    if (snesreason > -1) {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
		       "Solution time of %f sec, SNES converged %d (%s) \n",
		       v2 - v1, snesreason, SNESConvergedReasons[snesreason]);CHKERRQ(ierr);
    } else {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
		       "Solution time of %f sec, SNES diverged %d (%s).\n",
		       v2 - v1, snesreason, SNESConvergedReasons[snesreason]);CHKERRQ(ierr);
    }

    user->current_step++;

  }/* end while*/
  PetscFunctionReturn(0);
}

  std::pair<unsigned int, Real>
  PetscDMNonlinearSolver<T>::solve (SparseMatrix<T>& jac_in,  // System Jacobian Matrix
				    NumericVector<T>& x_in,   // Solution vector
				    NumericVector<T>& r_in,   // Residual vector
				    const double,             // Stopping tolerance
				    const unsigned int)
  {
    START_LOG("solve()", "PetscNonlinearSolver");
    this->init ();

    // Make sure the data passed in are really of Petsc types
    libmesh_cast_ptr<PetscMatrix<T>*>(&jac_in);
    libmesh_cast_ptr<PetscVector<T>*>(&r_in);

    // Extract solution vector
    PetscVector<T>* x = libmesh_cast_ptr<PetscVector<T>*>(&x_in);

    int ierr=0;
    int n_iterations =0;

    // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
    Real final_residual_norm=0.;

    if (this->user_presolve)
      this->user_presolve(this->system());

    //Set the preconditioning matrix
    if (this->_preconditioner)
      this->_preconditioner->set_matrix(jac_in);

    ierr = SNESSolve (this->_snes, PETSC_NULL, x->vec());
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

    ierr = SNESGetIterationNumber(this->_snes,&n_iterations);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

    ierr = SNESGetLinearSolveIterations(this->_snes, &this->_n_linear_iterations);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

    ierr = SNESGetFunctionNorm(this->_snes,&final_residual_norm);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

    // Get and store the reason for convergence
    SNESGetConvergedReason(this->_snes, &this->_reason);

    //Based on Petsc 2.3.3 documentation all diverged reasons are negative
    this->converged = (this->_reason >= 0);

    this->clear();

    STOP_LOG("solve()", "PetscNonlinearSolver");

    // return the # of its. and the final residual norm.
    return std::make_pair(n_iterations, final_residual_norm);
  }

int main(int argc, char **argv)
{
  SNES           snes;   /* nonlinear solver */
  DM             da;     /* grid */
  Vec            u;      /* solution vector */
  AppCtx         user;   /* user-defined work context */
  PetscReal      t;      /* time */
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, PETSC_NULL, help);CHKERRQ(ierr);
  /* Create solver */
  ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
  /* Create mesh */
  ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-4,3,1,PETSC_NULL,&da);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da, &user);CHKERRQ(ierr);
  ierr = SNESSetDM(snes, da);CHKERRQ(ierr);
  /* Create coefficient */
  ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-4,1,1,PETSC_NULL,&user.cda);CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(user.cda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(user.cda, &user.Kappa);CHKERRQ(ierr);
  ierr = FormPermeability(user.cda, user.Kappa, &user);CHKERRQ(ierr);
  /* Setup Problem */
  ierr = DMDASetLocalFunction(da, (DMDALocalFunction1) FormFunctionLocal);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(da, &u);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(da, &user.uold);CHKERRQ(ierr);
  user.sl  = 1.0;
  user.vl  = 0.1;
  user.pl  = 1.0;
  user.phi = 1.0;
  user.kappaWet   = 1.0;
  user.kappaNoWet = 0.3;
  /* Time Loop */
  user.dt = 0.1;
  for(PetscInt n = 0; n < 100; ++n, t += user.dt) {
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Starting time %g\n", t);CHKERRQ(ierr);
    ierr = VecView(u, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
    /* Solve */
    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
    ierr = SNESSolve(snes, PETSC_NULL, u);CHKERRQ(ierr);
    /* Update */
    ierr = VecCopy(u, user.uold);CHKERRQ(ierr);

    ierr = VecView(u, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
  }
  /* Cleanup */
  ierr = DMRestoreGlobalVector(da, &u);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(da, &user.uold);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(user.cda, &user.Kappa);CHKERRQ(ierr);
  ierr = DMDestroy(&user.cda);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

int main(int argc, char **argv)
{
  SNES           snes; /* nonlinear solver */
  DM             dm;   /* problem definition */
  Vec            u, r; /* solution and residual */
  AppCtx         user; /* user-defined work context */
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
  ierr = PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);CHKERRQ(ierr);
  ierr = SetupParameters(&user);CHKERRQ(ierr);
  ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
  ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr);
  /* Setup problem */
  ierr = PetscMalloc(2 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);CHKERRQ(ierr);
  ierr = SetupDiscretization(dm, &user);CHKERRQ(ierr);
  ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
  ierr = VecDuplicate(u, &r);CHKERRQ(ierr);

  ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr);

  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  {
    Parameter *param;
    void      *ctxs[2];

    ierr = PetscBagGetData(user.bag, (void **) &param);CHKERRQ(ierr);
    ctxs[0] = ctxs[1] = param;
    ierr = DMProjectFunction(dm, 0.0, user.exactFuncs, ctxs, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) u, "Exact Solution");CHKERRQ(ierr);
    ierr = VecViewFromOptions(u, NULL, "-exact_vec_view");CHKERRQ(ierr);
  }
  ierr = DMSNESCheckFromOptions(snes, u, NULL, NULL);CHKERRQ(ierr);
  ierr = VecSet(u, 0.0);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) u, "Solution");CHKERRQ(ierr);
  ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
  ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr);

  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&r);CHKERRQ(ierr);
  ierr = PetscFree(user.exactFuncs);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  ierr = PetscBagDestroy(&user.bag);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

/*
   Defines the RHS function that is passed to the time-integrator. 

   Solves F(U,V) for V and then computes f(U,V)

*/
PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
{
  AppCtx         *ctx = (AppCtx*)actx;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ctx->t = t;
  ctx->U = U;
  ierr   = SNESSolve(ctx->snes,NULL,ctx->V);CHKERRQ(ierr);
  ierr   = (*ctx->f)(t,U,ctx->V,F);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}