#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: gsnes.c,v 1.11 2000/07/16 05:40:27 knepley Exp $";
#endif

#include "src/snes/snesimpl.h"    /*I "petscsnes.h" I*/
#include "src/gvec/gvecimpl.h"    /*I "gsolver.h" I*/
#include "gsolver.h"

/*
  This file provides routines for grid SNES objects. They support solution of
  systems of nonlinear equations generated from a mesh and discretization
  specified by a Grid object.
*/

#undef  __FUNCT__
#define __FUNCT__ "GSNESCreateStructures_Private"
int GSNESCreateStructures_Private(GSNES snes, void *funcCtx, void *jacCtx)
{
  Grid       grid;
  GMat       J;
  GVec       f;
  PetscTruth isMatrixFree;
  int        ierr;

  PetscFunctionBegin;
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  ierr = GVecCreateConstrained(grid, &f);                                                                 CHKERRQ(ierr);
  ierr = GridIsMatrixFree(grid, &isMatrixFree);                                                           CHKERRQ(ierr);
  if (isMatrixFree == PETSC_TRUE) {
    if (grid->explicitConstraints == PETSC_TRUE) {
      ierr = GMatCreateMF(grid, grid->constraintOrder, grid->constraintOrder, &J);                        CHKERRQ(ierr);
    } else {
      ierr = GMatCreateMF(grid, grid->order, grid->order, &J);                                            CHKERRQ(ierr);
    }
  } else {
    if (grid->explicitConstraints == PETSC_TRUE) {
      ierr = GMatCreateRectangular(grid, grid->constraintOrder, grid->constraintOrder, &J);               CHKERRQ(ierr);
    } else {
      ierr = GMatCreate(grid, &J);                                                                        CHKERRQ(ierr);
    }
  }
  ierr = SNESSetFunction(snes, f, (int (*)(SNES, Vec, Vec, void *)) GSNESEvaluateRhs, funcCtx);           CHKERRQ(ierr);
  ierr = SNESSetJacobian(snes, J, J, (int (*)(SNES, Vec, Mat *, Mat *, MatStructure *, void *)) GSNESEvaluateJacobian, jacCtx);CHKERRQ(ierr);
  PetscLogObjectParent(grid, f);
  PetscLogObjectParent(grid, J);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESDestroyStructures_Private"
int GSNESDestroyStructures_Private(GSNES snes)
{
  GMat J;
  GVec f;
  int  ierr;

  PetscFunctionBegin;
  ierr = SNESGetFunction(snes, &f, PETSC_NULL, PETSC_NULL);                                               CHKERRQ(ierr);
  ierr = SNESGetJacobian(snes, &J, PETSC_NULL, PETSC_NULL, PETSC_NULL);                                   CHKERRQ(ierr);
  ierr = GMatDestroy(J);                                                                                  CHKERRQ(ierr);
  ierr = GVecDestroy(f);                                                                                  CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESDestroy"
/*@C
  GSNESDestroy - Destroys a grid SNES.

  Collective on GSNES

  Input Parameter:
. snes - The nonlinear solver context

  Level: beginner

.keywords: grid SNES, destroy
.seealso: SNESDestroy(), GSNESDuplicate()
@*/
int GSNESDestroy(GSNES snes)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  if (--snes->refct > 0)
    PetscFunctionReturn(0);
  ierr = GSNESDestroyStructures_Private(snes);                                                            CHKERRQ(ierr);
  ierr = SNESDestroy(snes);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESView"
/*@ 
  GSNESView - Views a grid SNES.

  Collective on GSNES

  Input Parameters:
+ snes   - The grid SNES
- viewer - [Optional] A visualization context

  Options Database Key:
$ -snes_view : calls SNESView() at end of SNESSolve()

  Notes:
  The available visualization contexts include
$   VIEWER_STDOUT_SELF - standard output (default)
$   VIEWER_STDOUT_WORLD - synchronized standard
$     output where only the first processor opens
$     the file.  All other processors send their 
$     data to the first processor to print. 

  The user can open alternative vistualization contexts with
$   PetscViewerFileOpenASCII() - output vector to a specified file

  Level: beginner

.keywords: grid SNES, view, visualize
.seealso: SNESView()
@*/
int GSNESView(GSNES snes, PetscViewer viewer)
{  
  Grid       grid;
  FILE      *fd;
  int        numActiveFields;
  int        field, func, op;
  PetscTruth isascii;
  int        ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  if (!viewer) {
    viewer = PETSC_VIEWER_STDOUT_SELF;
  } else {
    PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE);
  }

  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii);                            CHKERRQ(ierr);
  if (isascii == PETSC_TRUE) {
    ierr = GridView(grid, viewer);                                                                        CHKERRQ(ierr);
    ierr = PetscViewerFlush(viewer);                                                                      CHKERRQ(ierr);
    ierr = PetscViewerASCIIGetPointer(viewer, &fd);                                                       CHKERRQ(ierr);
    PetscFPrintf(grid->comm, fd, "GSNES Object:\n");
    ierr = GridGetNumActiveFields(grid, &numActiveFields);                                                CHKERRQ(ierr); 
    if (numActiveFields == 1) {
      PetscFPrintf(snes->comm, fd, "  %d active field:\n", numActiveFields);
    } else {
      PetscFPrintf(snes->comm, fd, "  %d active fields:\n", numActiveFields);
    }
    for(field = 0; field < grid->numFields; field++) {
      if (grid->fields[field].isActive == PETSC_FALSE) continue;
      if (grid->fields[field].name) {
        PetscFPrintf(grid->comm, fd, "  %s field with %d components:\n", grid->fields[field].name, grid->fields[field].numComp);
      } else {
        PetscFPrintf(grid->comm, fd, "  field %d with %d components:\n", field, grid->fields[field].numComp);
      }
      for(func = 0; func < grid->numRhsFuncs; func++) {
        if (grid->rhsFuncs[func].field == field) {
          PetscFPrintf(grid->comm, fd, "    Rhs function with scalar multiplier %g\n", PetscRealPart(grid->rhsFuncs[func].alpha));
        }
      }
      for(op = 0; op < grid->numRhsOps; op++) {
        if (grid->rhsOps[op].field == field) {
          if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) {
            PetscFPrintf(grid->comm, fd, "    Rhs nonlinear operator with scalar multiplier %g\n",
                         PetscRealPart(grid->rhsOps[field].alpha));
          } else {
            PetscFPrintf(grid->comm, fd, "    Rhs operator %d with scalar multiplier %g\n",
                         grid->rhsOps[op].op, PetscRealPart(grid->rhsOps[op].alpha));
          }
        }
      }
      for(op = 0; op < grid->numMatOps; op++) {
        if (grid->matOps[op].field == field) {
          PetscFPrintf(grid->comm, fd, "    Jacobian operator %d with scalar multiplier %g\n",
                       grid->matOps[op].op, PetscRealPart(grid->matOps[op].alpha));
        }
      }
    }
    ierr = SNESView(snes, viewer);                                                                        CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESDuplicate"
/*@C
  GSNESDuplicate - Duplicates a grid SNES.

  Collective on GSNES

  Input Parameter:
. snes    - The GSNES

  Output Parameter:
. newsnes - The new GSNES

  Level: beginner

.keywords: grid SNES, duplicate
.seealso: SNESDuplicate(), GSNESDestroy()
@*/
int GSNESDuplicate(GSNES snes, GSNES *newsnes)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  PetscValidPointer(newsnes);
  SETERRQ(PETSC_ERR_SUP, " ");
#if 0
  PetscFunctionReturn(0);
#endif
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESDuplicateMonitors"
/*@C
  GSNESDuplicateMonitors - Duplicates the monitors of a grid SNES in another.

  Collective on GSNES

  Input Parameter:
. snes    - The GSNES to copy from

  Output Parameter:
. newsnes - The GSNES to copy to

  Level: intermediate

.keywords: grid snes, duplicate, monitor
.seealso: GSNESDuplicate(), GSNESDestroy()
@*/
int GSNESDuplicateMonitors(GSNES snes, GSNES newsnes)
{
  int mon, mon2;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes,    SNES_COOKIE);
  PetscValidHeaderSpecific(newsnes, SNES_COOKIE);

  for(mon = 0; mon < snes->numbermonitors; mon++)
  {
    for(mon2 = 0; mon2 < newsnes->numbermonitors; mon2++)
      if (newsnes->monitor[mon] == snes->monitor[mon])
        break;
    if (mon2 == newsnes->numbermonitors) {
      ierr = SNESSetMonitor(newsnes, snes->monitor[mon], snes->monitorcontext[mon], PETSC_NULL);          CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESReallocate"
/*@C
  GSNESReallocate - Reallocates the storage for a grid SNES.

  Collective on GSNES

  Input Parameter:
. snes - The GSNES

  Level: advanced

.keywords: grid snes, reallocate
.seealso: GridReform(), GSNESDuplicate(), GSNESDestroy()
@*/
int GSNESReallocate(GSNES snes)
{
  void *funcCtx;
  void *jacCtx;
  int   ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  /* Preserve contexts */
  ierr = SNESGetFunction(snes, PETSC_NULL, &funcCtx, PETSC_NULL);                                         CHKERRQ(ierr);
  ierr = SNESGetJacobian(snes, PETSC_NULL, PETSC_NULL, &jacCtx, PETSC_NULL);                              CHKERRQ(ierr);
  ierr = GSNESDestroyStructures_Private(snes);                                                            CHKERRQ(ierr);
  ierr = GSNESCreateStructures_Private(snes, funcCtx, jacCtx);                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESGetGrid"
/*@
  GSNESGetGrid - Gets the grid from a grid SNES.

  Not collective

  Input Parameter:
. snes - The grid SNES

  Output Parameter:
. grid - The grid context

  Level: intermediate

.keywords: grid SNES, grid, get
.seealso: GVecGetGrid(), GMatGetGrid()
@*/
int GSNESGetGrid(GSNES snes, Grid *grid)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  PetscValidPointer(grid);

  ierr = PetscObjectQuery((PetscObject) snes, "Grid", (PetscObject *) grid);                              CHKERRQ(ierr);
  PetscValidHeaderSpecific(*grid, GRID_COOKIE);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateRhs"
/*@
  GSNESEvaluateRhs - Constructs the Rhs vector for Newton's equation.

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. x    - The current iterate
- ctx  - The optional user context

  Output Parameter:
. f    - The function value

  Note:
  This function initializes the vector f to zero before calculation.

  Level: advanced

.keywords: grid SNES, rhs
.seealso: GridEvaluateRhs()
@*/
int GSNESEvaluateRhs(GSNES snes, GVec x, GVec f, PetscObject ctx)
{
  Grid        grid;
  PetscScalar zero = 0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize vector */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Form function */
  ierr = GridEvaluateRhs(grid, x, f, ctx);                                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateRhsFunction"
/*@
  GSNESEvaluateRhsFunction - Constructs only the constant portion of the nonlinear equation,
  that is the portion arising from functions independent of the solution variable x.

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. x    - The current iterate
- ctx  - The optional user context

  Output Parameter:
. f    - The function value

  Notes:
  This function actually constructs -f since f is usually written on the right, and the user
  cannot just multiply by -1 after boundary conditions are imposed.

  Level: advanced

.keywords: grid SNES, rhs, function
.seealso: GSNESEvaluateRhsOperator(), GSNESEvaluateRhs(), GSNESEvaluateJacobian()
@*/
int GSNESEvaluateRhsFunction(GSNES snes, GVec x, GVec f, void *ctx)
{
  Grid        grid;
  PetscScalar zero = 0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize vector */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Form function */
  ierr = GridScaleRhs(grid, -1.0);                                                                        CHKERRQ(ierr);
  ierr = GridEvaluateRhsFunction(grid, x, f, ctx);                                                        CHKERRQ(ierr);
  ierr = GridScaleRhs(grid, -1.0);                                                                        CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateRhsOperator"
/*@
  GSNESEvaluateRhsOperator - Constructs only the portion of the nonlinear equation dependent of the solution variable x.

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. x    - The current iterate
- ctx  - The optional user context

  Output Parameter:
. f    - The function value

  Level: advanced

.keywords: grid SNES, rhs, operator
.seealso: GSNESEvaluateRhsLinearOperator(), GSNESEvaluateRhsNonlinearOperator(), GSNESEvaluateRhsFunction(),
          GSNESEvaluateRhs(), GSNESEvaluateJacobian()
@*/
int GSNESEvaluateRhsOperator(GSNES snes, GVec x, GVec f, void *ctx)
{
  Grid        grid;
  PetscTruth  reduceElement;
  PetscScalar zero = 0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize vector */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Turn off boundary values */
  ierr = GridGetReduceElement(grid, &reduceElement);                                                      CHKERRQ(ierr);
  ierr = GridSetReduceElement(grid, PETSC_FALSE);                                                         CHKERRQ(ierr);
  /* Apply operator */
  ierr = GridEvaluateRhsOperator(grid, x, f, PETSC_TRUE, PETSC_TRUE, ctx);                                CHKERRQ(ierr);
  /* Restore boundary values */
  ierr = GridSetReduceElement(grid, reduceElement);                                                       CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateRhsLinearOperator"
/*@
  GSNESEvaluateRhsLinearOperator - Constructs only the portion of the nonlinear equation
  linearly dependent on the solution variable x.

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. x    - The current iterate
- ctx  - The optional user context

  Output Parameter:
. f    - The function value

  Level: advanced

.keywords: grid SNES, rhs, operator, linear
.seealso: GSNESEvaluateRhsOperator(), GSNESEvaluateRhsFunction(), GSNESEvaluateRhs(), GSNESEvaluateJacobian()
@*/
int GSNESEvaluateRhsLinearOperator(GSNES snes, GVec x, GVec f, void *ctx)
{
  Grid        grid;
  PetscTruth  reduceElement;
  PetscScalar zero = 0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize vector */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Turn off boundary values */
  ierr = GridGetReduceElement(grid, &reduceElement);                                                      CHKERRQ(ierr);
  ierr = GridSetReduceElement(grid, PETSC_FALSE);                                                         CHKERRQ(ierr);
  /* Apply operator */
  ierr = GridEvaluateRhsOperator(grid, x, f, PETSC_TRUE, PETSC_FALSE, ctx);                               CHKERRQ(ierr);
  /* Restore boundary values */
  ierr = GridSetReduceElement(grid, reduceElement);                                                       CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateRhsNonlinearOperator"
/*@
  GSNESEvaluateRhsNonlinearOperator - Constructs only the portion of the nonlinear equation
  nonlinearly dependent on the solution variable x.

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. n    - The number of input vectors
. vecs - The input vectors
- ctx  - The optional user context

  Output Parameter:
. f    - The function value

  Level: advanced

  Note:
  If there is only one input argument, it will be used as every argument for
  multilinear functions.

.keywords: grid SNES, rhs, operator, nonlinear
.seealso: GSNESEvaluateRhsOperator(), GSNESEvaluateRhsFunction(), GSNESEvaluateRhs(), GSNESEvaluateJacobian()
@*/
int GSNESEvaluateRhsNonlinearOperator(GSNES snes, int n, GVec vecs[], GVec f, void *ctx)
{
  Grid              grid;
  NonlinearOperator op;
  int               field;
  PetscScalar       alpha;
  PetscTruth        isALE;
  PetscTruth        reduceElement;
  PetscScalar       zero = 0.0;
  int               ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize vector */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Turn off boundary values */
  ierr = GridGetReduceElement(grid, &reduceElement);                                                      CHKERRQ(ierr);
  ierr = GridSetReduceElement(grid, PETSC_FALSE);                                                         CHKERRQ(ierr);
  /* Apply operator */
  if (n == 1) {
    ierr = GridEvaluateRhsOperator(grid, vecs[0], f, PETSC_FALSE, PETSC_TRUE, ctx);                       CHKERRQ(ierr);
  } else {
    ierr = GridGetNonlinearOperatorStart(grid, &op, &field, &alpha, &isALE);                              CHKERRQ(ierr);
    while(op != PETSC_NULL) {
      ierr = GVecEvaluateNonlinearOperatorGalerkin(f, n, vecs, 1, &field, op, alpha, isALE, ctx);         CHKERRQ(ierr);
      ierr = GridGetNonlinearOperatorNext(grid, &op, &field, &alpha, &isALE);                             CHKERRQ(ierr);
    }
  }
  /* Restore boundary values */
  ierr = GridSetReduceElement(grid, reduceElement);                                                       CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateJacobian"
/*@
  GSNESEvaluateJacobian - Constructs the system matrix for Newton's equation

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. x    - The current iterate
. flag - The matrix structure flag
- ctx  - The optional user context

  Output Parameters:
+ J    - The Jacobian matrix
- M    - The preconditioner for the Jacobian matrix, usually the same as J

  Level: advanced

.keywords: grid SNES, jacobian
.seealso: GSNESEvaluateJacobianMF(), GSNESEvaluateRhs(), GridEvaluateSystemMatrix()
@*/
int GSNESEvaluateJacobian(GSNES snes, GVec x, GMat *J, GMat *M, MatStructure *flag, PetscObject ctx)
{
  Grid grid;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize matrix */
  ierr = MatZeroEntries(*J);                                                                              CHKERRQ(ierr);
  /* Form function */
  ierr = GridEvaluateSystemMatrix(grid, x, J, M, flag, ctx);                                              CHKERRQ(ierr);
  /* Apply boundary conditions */
  ierr = GMatSetBoundary(*J, 1.0, ctx);                                                                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESEvaluateJacobianMF"
/*@
  GSNESEvaluateJacobianMF - Constructs the application of the system matrix for Newton's equation

  Collective on GSNES

  Input Parameters:
+ snes - The grid SNES
. x    - The current iterate
. y    - The vector to which to apply the Jacobian
- ctx  - The optional user context

  Output Parameter:
. f    - The vector J(x) y

  Level: advanced

.keywords: grid SNES< jacobian, mf
.seealso: GSNESEvaluateJacobian(), GSNESEvaluateRhs(), GridEvaluateSystemMatrix(), GVecEvaluateJacobian()
@*/
int GSNESEvaluateJacobianMF(GSNES snes, GVec x, GVec y, GVec f, void *ctx)
{
  Grid        grid;
  PetscScalar zero = 0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  /* Initialize vector  */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Form function */
  ierr = GVecEvaluateJacobian(f, x, y, ctx);                                                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GSNESRhsBC"
/*@
  GSNESRhsBC - This sets zero boundary conditions for the Rhs.

  Collective on GSNES

  Input Parameters:
+ snes - The SNES
- ctx  - The user-supplied context

  Output Parameter:
. rhs  - The modified Rhs

  Level: advanced

.keywords: grid SNES, set, boundary conditions
.seealso GSNESSetSolutionBC()
@*/
int GSNESRhsBC(GSNES snes, GVec rhs, void *ctx)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  PetscValidHeaderSpecific(rhs,  VEC_COOKIE);
  ierr = GVecSetBoundaryZero(rhs, ctx);                                                                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GSNESSolutionBC"
/*@
  GSNESSolutionBC - This sets the boundary conditions registered with the grid
  on the solution vector.

  Collective on GSNES

  Input Parameters:
+ snes - The SNES
- ctx  - The user-supplied context

  Output Parameter:
. sol  - The modified solution

  Level: advanced

.keywords: GSNES, set, boundary conditions
.seealso: GSNESSetRhsBC()
@*/
int GSNESSolutionBC(GSNES snes, GVec sol, void *ctx)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  PetscValidHeaderSpecific(sol,  VEC_COOKIE);
  ierr = GVecSetBoundary(sol, ctx);                                                                       CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GSNESUpdate"
/*@
  GSNESUpdate - This is a general purpose function which is called at the beginning of every iterate.

  Collective on GSNES

  Input Parameters:
+ snes - The nonlinear solver context
- step - The current step

  Level: advanced

.keywords: GSNES, update
.seealso: SNESSolve()
@*/
int GSNESUpdate(GSNES snes, int step)
{
  Grid       grid;
  Mesh       mesh;
  MeshMover  mover;
  PetscTruth isMoving;
  int        ierr;

  PetscFunctionBegin;
  /* Calculate mesh velocity */
  ierr = GSNESGetGrid(snes, &grid);                                                                       CHKERRQ(ierr);
  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  ierr = MeshGetMovement(mesh, &isMoving);                                                                CHKERRQ(ierr);
  if (isMoving == PETSC_TRUE) {
    ierr = MeshGetMover(mesh, &mover);                                                                    CHKERRQ(ierr);
    ierr = MeshMoverCalcNodeVelocities(mover, PETSC_TRUE);                                                CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESCreate"
/*@
  GSNESCreate - Creates a grid SNES associated with a particular discretization context.

  Collective on Grid

  Input Parameter:
+ grid - The discretization context
- ctx  - The optional user context

  Output Parameter:
. snes - The nonlinear solver context

  Level: beginner

.keywords: GSNES, create
.seealso: GSNESDestroy()
@*/
int GSNESCreate(Grid grid, void *ctx, GSNES *snes)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(snes);
#ifndef PETSC_USE_DYNAMIC_LIBRARIES
  ierr = GSolverInitializePackage(PETSC_NULL);                                                            CHKERRQ(ierr);
#endif

  /* Create the SNES */
  ierr = SNESCreate(grid->comm,snes);                                                                     CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject) *snes, "Grid", (PetscObject) grid);                             CHKERRQ(ierr);
  (*snes)->isGSNES = PETSC_TRUE;

  /* Set boundary conditions */
  ierr = SNESSetRhsBC(*snes, GSNESRhsBC);                                                                 CHKERRQ(ierr);
  ierr = SNESSetSolutionBC(*snes, GSNESSolutionBC);                                                       CHKERRQ(ierr);
  if (grid->reduceSystem == PETSC_TRUE) {
    /* Save space and time by reducing the system right in construction rather than using an explicit matvec and subtraction */
    grid->reduceElement = PETSC_TRUE;
    /* Since we solve J x = -F, we need a BC multiplier of -1 */
    ierr = GridSetBCMultiplier(grid, -1.0);                                                               CHKERRQ(ierr);
  }

  /* Set update function */
  ierr = SNESSetUpdate(*snes, GSNESUpdate);                                                               CHKERRQ(ierr);

  /* Setup SNES data structures */
  ierr = GSNESCreateStructures_Private(*snes, ctx, ctx);                                                  CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESDefaultComputeJacobianWithColoring"
int GSNESDefaultComputeJacobianWithColoring(SNES snes, GVec x, GMat *J, GMat *M, MatStructure *flag, void *ctx)
{
  SETERRQ(PETSC_ERR_SUP, " ");
#if 0
  Grid grid;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  PetscValidHeaderSpecific(x,    VEC_COOKIE);
  PetscValidHeaderSpecific(J,    MAT_COOKIE);
  PetscValidHeaderSpecific(M,    MAT_COOKIE);
  ierr = GMatGetGrid(*M, &grid);                                                                          CHKERRQ(ierr);
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  if (!grid->fdcoloring) {
    ierr = GMatFDColoringCreate(*M);                                                                      CHKERRQ(ierr);
  }
  ierr = SNESDefaultComputeJacobianWithColoring(snes, x, J, M, flag, grid->fdcoloring);                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
#endif
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESSolutionMonitor"
/*@
  GSNESSolutionMonitor - Monitors solution at each GSNES iteration.

  Collective on GSNES

  Input Parameters:
+ snes  - The nonlinear solver context
. it    - The number of iterations so far
. rnorm - The current (approximate) residual norm
- ctx   - The viewer

  Level: advanced

.keywords: grid SNES, monitor, solution
.seealso: GSNESResidualMonitor(), GSNESErrorMonitor(), SNESDefaultMonitor(), SNESSetMonitor()
@*/
int GSNESSolutionMonitor(GSNES snes, int it, PetscReal rnorm, void *ctx)
{
  PetscViewer viewer = (PetscViewer) ctx;
  Vec         x;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = SNESGetSolution(snes, &x);                                                                       CHKERRQ(ierr);
  ierr = VecView(x, viewer);                                                                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESResidualMonitor"
/*@
  GSNESResidualMonitor - Monitors residual at each SNES iteration.

  Collective on GSNES

  Input Parameters:
+ snes  - The nonlinear solver context
. it    - The number of iterations so far
. rnorm - The current (approximate) residual norm
- ctx   - The viewer

  Level: advanced

.keywords: grid SNES, monitor, residual
.seealso: GSNESSolutionMonitor(), GSNESErrorMonitor(), SNESDefaultMonitor(), SNESSetMonitor()
@*/
int GSNESResidualMonitor(GSNES snes, int it, PetscReal rnorm, void *ctx)
{
  PetscViewer viewer = (PetscViewer) ctx;
  Vec         x;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  ierr = SNESGetFunction(snes, &x, PETSC_NULL, PETSC_NULL);                                               CHKERRQ(ierr);
  ierr = VecView(x, viewer);                                                                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GSNESErrorMonitor"
/*@
  GSNESErrorMonitor - Displays the error at each iteration.

  Collective on GSNES

  Input Parameters:
+ snes   - The nonlinear context
. it     - The number of iterations so far
. rnorm  - The current (approximate) residual norm
- monCtx - A GSNESErrorMonitorCtx

  Notes: 
  The final argument to SNESSetMonitor() with this routine must be a pointer to a GSNESErrorMonitorCtx.

  Level: advanced

.keywords: grid SNES, monitor, residual
.seealso: GSNESSolutionMonitor(), GSNESResidualMonitor(), SNESDefaultMonitor(), SNESSetMonitor()
@*/
int GSNESErrorMonitor(GSNES snes, int it, PetscReal rnorm, void *monCtx)
{
  GSNESErrorMonitorCtx *errorCtx = (GSNESErrorMonitorCtx *) monCtx;
  void                 *ctx      = errorCtx->ctx;
  PetscScalar           minusOne = -1.0;
  GVec                  x, e;
  PetscReal             norm_2, norm_max, sol_norm_2;
  int                   ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_COOKIE);
  if (it == 0) {
    /* Build solution at the beginning */
    ierr = (*errorCtx->solutionFunc)(errorCtx->solution, ctx);                                            CHKERRQ(ierr);
  }
  ierr = SNESGetSolution(snes, &x);                                                                       CHKERRQ(ierr);
  ierr = GVecGetWorkGVec(errorCtx->solution, &e);                                                         CHKERRQ(ierr);
  ierr = VecWAXPY(&minusOne, x, errorCtx->solution, e);                                                   CHKERRQ(ierr);
  ierr = VecView(e, errorCtx->error_viewer);                                                              CHKERRQ(ierr);

  /* Compute 2-norm and max-norm of error */
  if (errorCtx->norm_error_viewer) {
    ierr = VecNorm(e, NORM_2, &norm_2);                                                                   CHKERRQ(ierr);
    ierr = VecNorm(e, NORM_MAX, &norm_max);                                                               CHKERRQ(ierr);
    ierr = VecNorm(errorCtx->solution, NORM_2, &sol_norm_2);                                              CHKERRQ(ierr);
    PetscViewerASCIIPrintf(errorCtx->norm_error_viewer,
                           "Iteration %d residual norm %g error 2-norm %g error max-norm %g exact sol norm-2 %g\n",
                           it, rnorm, norm_2, norm_max, sol_norm_2);
  }
  ierr = GVecRestoreWorkGVec(errorCtx->solution, &e);                                                     CHKERRQ(ierr);  
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "GSNESOptionsChecker_Private"
int GSNESOptionsChecker_Private(GSNES snes) 
{
  char       *prefix;  
  int         loc[4], nmax;
  MPI_Comm    comm;
  PetscViewer viewer;
  PetscDraw   draw;
  PetscTruth  opt;
  int         ierr;

  PetscFunctionBegin;
  ierr = SNESGetOptionsPrefix(snes, &prefix);                                                             CHKERRQ(ierr);
  ierr = PetscObjectGetComm((PetscObject) snes, &comm);                                                   CHKERRQ(ierr);

  nmax = 4;
  loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
  ierr = PetscOptionsGetIntArray(prefix, "-gvec_snes_solutionmonitor", loc, &nmax, &opt);                 CHKERRQ(ierr);
  if (opt) {
    ierr = PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);                      CHKERRQ(ierr);
    ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);                                                      CHKERRQ(ierr);
    ierr = PetscDrawSetTitle(draw, "Approx. Solution");                                                   CHKERRQ(ierr);
    PetscLogObjectParent(snes, (PetscObject) viewer);
    ierr = SNESSetMonitor(snes, GSNESSolutionMonitor, (void *) viewer, PETSC_NULL);                       CHKERRQ(ierr);
  }
  nmax = 4;
  loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
  ierr = PetscOptionsGetIntArray(prefix, "-gvec_snes_residualmonitor", loc, &nmax, &opt);                 CHKERRQ(ierr);
  if (opt) {
    ierr = PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);                      CHKERRQ(ierr);
    ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);                                                      CHKERRQ(ierr);
    ierr = PetscDrawSetTitle(draw, "Residual");                                                           CHKERRQ(ierr);
    PetscLogObjectParent(snes, (PetscObject) viewer);
    ierr = SNESSetMonitor(snes, GSNESResidualMonitor, (void *) viewer, PETSC_NULL);                       CHKERRQ(ierr);
  }
#if 0
  nmax = 4;
  loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
  ierr = PetscOptionsGetIntArray(prefix, "-gvec_snes_errormonitor", loc, &nmax, &opt);                    CHKERRQ(ierr);
  if (opt) {
    ierr = PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);                      CHKERRQ(ierr);
    ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);                                                      CHKERRQ(ierr);
    ierr = PetscDrawSetTitle(draw, "Error");                                                              CHKERRQ(ierr);
    PetscLogObjectParent(snes, (PetscObject) viewer);
    ierr = SNESSetMonitor(snes, GSNESErrorMonitor, (void *) viewer, PETSC_NULL);                          CHKERRQ(ierr);
  }
#endif
  /*--------------------------------------------------------------------- */
  ierr = PetscOptionsHasName(prefix, "-gvec_snes_fd", &opt);                                              CHKERRQ(ierr);
  if (opt) {
    GMat A,B;
    ierr = SNESGetJacobian(snes, &A, &B, PETSC_NULL, PETSC_NULL);                                         CHKERRQ(ierr);
    ierr = SNESSetJacobian(snes, A, B, GSNESDefaultComputeJacobianWithColoring, 0);                       CHKERRQ(ierr);
    PetscLogInfo(snes, "GSNESSetFromOptions: Setting gvec finite difference Jacobian matrix\n");
  }
  /*--------------------------------------------------------------------- */
  ierr = PetscOptionsHasName(PETSC_NULL, "-help", &opt);                                                  CHKERRQ(ierr);
  if (opt) {
    char *pprefix;
    int   len = 2;
    if (prefix != PETSC_NULL) {
      ierr = PetscStrlen(prefix, &len);                                                                   CHKERRQ(ierr);
    }
    ierr = PetscMalloc((len+2) * sizeof(char), &pprefix);                                                 CHKERRQ(ierr);
    PetscStrcpy(pprefix, "-");
    if (prefix != PETSC_NULL)
      PetscStrcat(pprefix, prefix);
    PetscPrintf(comm," Additional SNES Monitor options for grid vectors\n");
    PetscPrintf(comm,"   %sgvec_snes_solutionmonitor\n", pprefix);
    PetscPrintf(comm,"   %sgvec_snes_residualmonitor\n", pprefix);
    PetscPrintf(comm,"   %sgvec_snes_errormonitor -- not working yet\n", pprefix);
    ierr = PetscFree(pprefix);                                                                            CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}
EXTERN_C_END