FiniteDifferencePreconditioner: now allow users to

specify finite difference type as coloring or standard.

If a standard FDP is used, matrix is set to full because
standard FDP will add off-diagonal entries

Refs #5819

framework/include/base/NonlinearSystem.h

@@ -73,6 +73,26 @@ class NonlinearSystem : public NonlinearSystemBase
 
 protected:
   TransientNonlinearImplicitSystem & _transient_sys;
+
+private:
+  /**
+  * Form preconditioning matrix via a standard finite difference method
+  * column-by-column. This method computes both diagonal and off-diagonal
+  * entrices regardless of the structure pattern of the Jacobian matrix.
+  */
+  void setupStandardFiniteDifferencedPreconditioner();
+
+  /**
+  * According to the nonzero pattern provided in the matrix, a graph is constructed.
+  * A coloring algorithm is applied to the graph. The graph is partitioned into several
+  * independent subgraphs (colors), and a finte difference method is applied color-by-color
+  * to form a preconditioning matrix. If the number of colors is small, this method is much
+  * faster than the standard FD. But there is an issue. If the matrix provided by users does not
+  * represent the actual structure of the true Jacobian, the matrix computed via coloring could
+  * be wrong or inaccurate. In this case, users should switch to the standard finite difference
+  * method.
+  */
+  void setupColoringFiniteDifferencedPreconditioner();
 };
 
 #endif /* NONLINEARSYSTEM_H */

framework/include/preconditioners/FiniteDifferencePreconditioner.h

@@ -16,6 +16,7 @@
 #define FINITEDIFFERENCEPRECONDITIONER_H
 
 #include "MoosePreconditioner.h"
+#include "MooseEnum.h"
 
 class FiniteDifferencePreconditioner;
 
@@ -29,6 +30,10 @@ class FiniteDifferencePreconditioner : public MoosePreconditioner
 {
 public:
   FiniteDifferencePreconditioner(const InputParameters & params);
+  MooseEnum & finiteDifferenceType() { return _finite_difference_type; }
+
+private:
+  MooseEnum _finite_difference_type;
 };
 
 #endif /* FINITEDIFFERENCEPRECONDITIONER_H */

framework/src/base/NonlinearSystem.C

@@ -16,6 +16,7 @@
 #include "NonlinearSystem.h"
 #include "FEProblem.h"
 #include "TimeIntegrator.h"
+#include "FiniteDifferencePreconditioner.h"
 
 #include "libmesh/nonlinear_solver.h"
 #include "libmesh/petsc_nonlinear_solver.h"
@@ -212,6 +213,44 @@ NonlinearSystem::stopSolve()
 
 void
 NonlinearSystem::setupFiniteDifferencedPreconditioner()
+{
+  std::shared_ptr<FiniteDifferencePreconditioner> fdp =
+      std::dynamic_pointer_cast<FiniteDifferencePreconditioner>(_preconditioner);
+  if (!fdp)
+    mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
+               "block with type = fdp");
+
+  if (fdp->finiteDifferenceType() == "coloring")
+    setupColoringFiniteDifferencedPreconditioner();
+  else if (fdp->finiteDifferenceType() == "standard")
+    setupStandardFiniteDifferencedPreconditioner();
+  else
+    mooseError("Unknown finite difference type");
+}
+
+void
+NonlinearSystem::setupStandardFiniteDifferencedPreconditioner()
+{
+#if LIBMESH_HAVE_PETSC
+  PetscNonlinearSolver<Number> * petsc_nonlinear_solver =
+      static_cast<PetscNonlinearSolver<Number> *>(_transient_sys.nonlinear_solver.get());
+
+  PetscMatrix<Number> * petsc_mat = static_cast<PetscMatrix<Number> *>(_transient_sys.matrix);
+
+  SNESSetJacobian(petsc_nonlinear_solver->snes(),
+                  petsc_mat->mat(),
+                  petsc_mat->mat(),
+#if PETSC_VERSION_LESS_THAN(3, 4, 0)
+                  SNESDefaultComputeJacobian,
+#else
+                  SNESComputeJacobianDefault,
+#endif
+                  nullptr);
+#endif
+}
+
+void
+NonlinearSystem::setupColoringFiniteDifferencedPreconditioner()
 {
 #ifdef LIBMESH_HAVE_PETSC
   // Make sure that libMesh isn't going to override our preconditioner

framework/src/preconditioners/FiniteDifferencePreconditioner.C

@@ -46,11 +46,18 @@ validParams<FiniteDifferencePreconditioner>()
                         "matrix for degrees of freedom that might be coupled "
                         "by inspection of the geometric search objects.");
 
+  MooseEnum finite_difference_type("standard coloring", "coloring");
+  params.addParam<MooseEnum>("finite_difference_type",
+                             finite_difference_type,
+                             "standard: standard finite difference"
+                             "coloring: finite difference based on coloring");
+
   return params;
 }
 
 FiniteDifferencePreconditioner::FiniteDifferencePreconditioner(const InputParameters & params)
-  : MoosePreconditioner(params)
+  : MoosePreconditioner(params),
+    _finite_difference_type(getParam<MooseEnum>("finite_difference_type"))
 {
   if (n_processors() > 1)
     mooseError("Can't use the Finite Difference Preconditioner in parallel yet!");
@@ -62,6 +69,10 @@ FiniteDifferencePreconditioner::FiniteDifferencePreconditioner(const InputParame
 
   bool full = getParam<bool>("full");
 
+  // standard finite difference method will add off-diagonal entries
+  if (_finite_difference_type == "standard")
+    full = true;
+
   if (!full)
   {
     // put 1s on diagonal