Merge branch 'fusion-ldrd-dev' into epic-dev

CMakeLists.txt

@@ -267,6 +267,15 @@ if (MFEM_USE_SUNDIALS)
   find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
 endif()
 
+# EPIC
+if (MFEM_USE_EPIC)
+  if (NOT (MFEM_USE_MPI AND MFEM_USE_SUNDIALS AND MFEM_USE_LAPACK) )
+    message(FATAL_ERROR " *** EPIC requires that MPI, SUNDIALS and LAPACK be enabled.")
+  else()
+    find_package(EPIC REQUIRED SUNDIALS NVector_Serial NVector_Parallel BLAS LAPACK)
+  endif()
+endif()
+
 # Mesquite
 if (MFEM_USE_MESQUITE)
   find_package(Mesquite REQUIRED)
@@ -428,7 +437,7 @@ endif()
 #    With newer versions of SuiteSparse which include METIS header using 64-bit
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
-set(MFEM_TPLS MPI_CXX OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS PETSC
+set(MFEM_TPLS MPI_CXX OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS EPIC PETSC
     SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB NETCDF
     MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2
     CUSPARSE MKL_CPARDISO AMGX CALIPER)

config/cmake/MFEMConfig.cmake.in

@@ -29,6 +29,7 @@ set(MFEM_USE_LEGACY_OPENMP @MFEM_USE_LEGACY_OPENMP@)
 set(MFEM_USE_MEMALLOC @MFEM_USE_MEMALLOC@)
 set(MFEM_TIMER_TYPE @MFEM_TIMER_TYPE@)
 set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
+set(MFEM_USE_EPIC @MFEM_USE_EPIC@)
 set(MFEM_USE_MESQUITE @MFEM_USE_MESQUITE@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)

config/cmake/config.hpp.in

@@ -165,6 +165,9 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 #cmakedefine MFEM_USE_SUNDIALS
 
+// Enable MFEM functionality based on the EPIC libraries.
+#cmakedefine MFEM_USE_EPIC
+
 // Version of HYPRE used for building MFEM.
 #cmakedefine MFEM_HYPRE_VERSION @MFEM_HYPRE_VERSION@
 

config/cmake/modules/FindEPIC.cmake

@@ -0,0 +1,21 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - EPIC_FOUND
+#   - EPIC_LIBRARIES
+#   - EPIC_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(EPIC EPIC EPIC_DIR
+       "include" Epic.h "lib" epic1.0.0
+  "Paths to headers required by EPIC." "Libraries required by EPIC.")
+

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -759,7 +759,7 @@ function(mfem_export_mk_files)
   set(CONFIG_MK_BOOL_VARS MFEM_USE_MPI MFEM_USE_METIS MFEM_USE_METIS_5
       MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB MFEM_USE_LIBUNWIND
       MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_OPENMP MFEM_USE_LEGACY_OPENMP
-      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
+      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_EPIC MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
       MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX
       MFEM_USE_GNUTLS MFEM_USE_GSLIB MFEM_USE_NETCDF MFEM_USE_PETSC
       MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT MFEM_USE_PUMI

config/config.hpp.in

@@ -85,6 +85,9 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 // #define MFEM_USE_SUNDIALS
 
+// Enable MFEM functionality based on the EPIC libraries.
+// #define MFEM_USE_EPIC
+
 // Enable MFEM functionality based on the Mesquite library.
 // #define MFEM_USE_MESQUITE
 

config/config.mk.in

@@ -29,6 +29,7 @@ MFEM_USE_OPENMP        = @MFEM_USE_OPENMP@
 MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
+MFEM_USE_EPIC          = @MFEM_USE_EPIC@
 MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@

config/defaults.cmake

@@ -30,6 +30,7 @@ option(MFEM_USE_OPENMP "Enable the OpenMP backend" OFF)
 option(MFEM_USE_LEGACY_OPENMP "Enable legacy OpenMP usage" OFF)
 option(MFEM_USE_MEMALLOC "Enable the internal MEMALLOC option." ON)
 option(MFEM_USE_SUNDIALS "Enable SUNDIALS usage" OFF)
+option(MFEM_USE_EPIC "Enable EPIC usage" OFF)
 option(MFEM_USE_MESQUITE "Enable MESQUITE usage" OFF)
 option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
@@ -115,6 +116,9 @@ set(SUNDIALS_DIR "${MFEM_DIR}/../sundials-5.0.0/instdir" CACHE PATH
 # set(SUNDIALS_REQUIRED_PACKAGES "SuiteSparse/KLU/AMD/BTF/COLAMD/config"
 #     CACHE STRING "Additional packages required by SUNDIALS.")
 
+set(EPIC_DIR "${MFEM_DIR}/../epic-cpp/instdir" CACHE PATH
+    "Path to the EPIC library.")
+
 set(MESQUITE_DIR "${MFEM_DIR}/../mesquite-2.99" CACHE PATH
     "Path to the Mesquite library.")
 

config/defaults.mk

@@ -122,6 +122,7 @@ MFEM_USE_LEGACY_OPENMP = NO
 MFEM_USE_MEMALLOC      = YES
 MFEM_TIMER_TYPE        = $(if $(NOTMAC),2,4)
 MFEM_USE_SUNDIALS      = NO
+MFEM_USE_EPIC          = NO
 MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
@@ -231,6 +232,11 @@ endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES
 
+# EPIC library configuration
+MESQUITE_DIR = @MFEM_DIR@/../epic-cpp/instdir
+MESQUITE_OPT = -I$(EPIC_DIR)/include
+MESQUITE_LIB = -L$(EPIC_DIR)/lib -lepic1.0.0
+
 # MESQUITE library configuration
 MESQUITE_DIR = @MFEM_DIR@/../mesquite-2.99
 MESQUITE_OPT = -I$(MESQUITE_DIR)/include

examples/CMakeLists.txt

@@ -159,6 +159,11 @@ if (MFEM_USE_AMGX)
   add_subdirectory(amgx)
 endif()
 
+# Include the examples/epic directory if EPIC is enabled.
+if (MFEM_USE_EPIC)
+  add_subdirectory(epic)
+endif()
+
 # Include the examples/ginkgo directory if GINKGO is enabled.
 if (MFEM_USE_GINKGO)
   add_subdirectory(ginkgo)

examples/epic/CMakeLists.txt

@@ -0,0 +1,64 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+set(EPIC_EXAMPLES_SRCS)
+list(APPEND EPIC_EXAMPLES_SRCS
+  ex16.cpp
+  )
+
+if (MFEM_USE_MPI)
+list(APPEND EPIC_EXAMPLES_SRCS
+    ex16p.cpp
+    )
+endif()
+
+# Include the source directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
+
+# Add "test_epic" target, see below.
+add_custom_target(test_epic
+  ${CMAKE_CTEST_COMMAND} -R epic USES_TERMINAL)
+
+# Add one executable per cpp file, adding "epic_" as prefix. Sets
+# "test_epic" as a target that depends on the given examples.
+set(PFX epic_)
+add_mfem_examples(EPIC_EXAMPLES_SRCS ${PFX} "" test_epic)
+
+# Testing.
+# The EPIC tests can be run separately using the target "test_epic"
+# which builds the examples and runs:
+#   ctest -R epic
+
+# Example 16: use the default options
+
+# Add the tests: one test per source file.
+foreach(SRC_FILE ${EPIC_EXAMPLES_SRCS})
+  get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
+  string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+  string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
+  set(TEST_NAME ${PFX}${TEST_NAME})
+
+  set(THIS_TEST_OPTIONS "-no-vis")
+  list(APPEND THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
+  # message(STATUS "Test ${TEST_NAME} options: ${THIS_TEST_OPTIONS}")
+
+  if (NOT (${TEST_NAME} MATCHES ".*p$"))
+    add_test(NAME ${TEST_NAME}_ser
+      COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
+  else()
+    add_test(NAME ${TEST_NAME}_np=4
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
+      ${MPIEXEC_POSTFLAGS})
+  endif()
+endforeach()
+

examples/epic/README

@@ -0,0 +1,17 @@
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               http://mfem.org
+
+This directory contains modifications of the example codes that illustrate the
+use of MFEM features based on the EPIC suite of time integration.
+
+To build these examples, make sure that MFEM is configured with the option
+"MFEM_USE_EPIC = YES".
+
+We recommend comparing the original example codes with the corresponding files
+in the current directory.

examples/epic/ex16.cpp

@@ -0,0 +1,610 @@
+//                                MFEM Example 16
+//                               EPIC Modification
+//
+// Compile with: make ex16
+//
+// Sample runs:  ex16
+//               ex16 -m ../../data/inline-tri.mesh
+//               ex16 -m ../../data/disc-nurbs.mesh -tf 2
+//               ex16 -s 8 -a 1.0 -k 0.0 -dt 1e-4 -tf 5e-2 -vs 25
+//               ex16 -m ../../data/fichera-q2.mesh
+//               ex16 -m ../../data/escher.mesh
+//               ex16 -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
+//               ex16 -m ../../data/amr-quad.mesh -o 4 -r 0
+//               ex16 -m ../../data/amr-hex.mesh -o 2 -r 0
+//
+// Description:  This example solves a time dependent nonlinear heat equation
+//               problem of the form du/dt = C(u), with a non-linear diffusion
+//               operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
+//
+//               We recommend viewing examples 2, 9 and 10 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ImplicitSolveOperator;
+class JacobianOperator;
+
+/** After spatial discretization, the conduction model can be written as:
+ *
+ *     du/dt = M^{-1}(-K(u) u)
+ *
+ *  where u is the vector representing the temperature, M is the mass matrix,
+ *  and K is the diffusion operator with diffusivity depending on u:
+ *  (\kappa + \alpha u).
+ *
+ *  Class ConductionOperator represents the right-hand side of the above ODE.
+ */
+class ConductionOperator : public TimeDependentOperator
+{
+protected:
+   FiniteElementSpace &fespace;
+   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
+
+   BilinearForm *M;
+   mutable BilinearForm *K;
+   mutable BilinearForm *dK;
+   mutable BilinearForm *J_K;
+
+   SparseMatrix Mmat;
+   mutable SparseMatrix J_K_mat;
+
+   mutable CGSolver M_solver; // Krylov solver for inverting the mass matrix M
+   DSmoother M_prec;  // Preconditioner for the mass matrix M
+
+   CGSolver Jg_solver; // Krylov solver for inverting the Jacobian in the nonlinear solve
+   DSmoother Jg_prec;  // Preconditioner for the Jacobian Jg
+
+   NewtonSolver newton_solver;
+   mutable JacobianOperator *jac;
+
+   double alpha, kappa;
+
+   mutable Vector z; // auxiliary vector
+
+   mutable int nRhsMult, nSetJac, nJacMult, nImpSolve, nImpIter, nImpMult, nImpSet;
+
+public:
+   Vector u0;
+
+   ConductionOperator(FiniteElementSpace &f, double alpha, double kappa, const Vector &u);
+
+   void UpdateStats();
+   void PrintStats(ostream& out);
+
+   void ExtractJacobians(const Vector& x, std::ostream &out, std::ostream &out2);
+
+   BilinearForm& GetKLambda(const Vector& u) const;
+   BilinearForm& GetdKLambda(const Vector& u) const;
+
+   virtual void Mult(const Vector &u, Vector &du_dt) const;
+   virtual Operator& GetGradient(const Vector &k) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
+
+   virtual ~ConductionOperator();
+};
+
+class ImplicitSolveOperator : public Operator
+{
+private:
+   double dt;
+   const Vector* x;
+   ConductionOperator* oper;
+
+   const SparseMatrix* M;
+   mutable SparseMatrix* Jg;
+
+   mutable Vector u, z;
+   mutable int nMult, nSet;
+
+public:
+   ImplicitSolveOperator(ConductionOperator* oper, const SparseMatrix* M, double dt, const Vector* x);
+
+   int GetnMult() { return nMult; }
+   int GetnSet() { return nSet; }
+   virtual void Mult(const Vector &k, Vector &gk) const;
+   virtual Operator &GetGradient(const Vector &k) const;
+};
+
+class JacobianOperator : public Operator
+{
+private:
+   Operator* J;
+   Operator* M_solver;
+
+   mutable int nMult;
+   mutable Vector z;
+public:
+   JacobianOperator(Operator* J, Operator* M_solver);
+
+   int GetnMult() { return nMult; }
+
+   void ExtractJacobian(const Vector& x, std::ostream &out);
+   virtual void Mult(const Vector &k, Vector &gk) const;
+};
+
+double InitialTemperature(const Vector &x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int ref_levels = 2;
+   int order = 2;
+   int ode_solver_type = 8; // Exponential Euler
+   double t_final = 0.5;
+   double dt = 1.0e-2;
+   double alpha = 1.0e-2;
+   double kappa = 0.5;
+   bool visualization = true;
+   bool visit = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver:\n\t"
+                  "1  - Forward Euler,\n\t"
+                  "2  - RK2,\n\t"
+                  "3  - RK3 SSP,\n\t"
+                  "4  - RK4,\n\t"
+                  "5  - Backward Euler,\n\t"
+                  "6  - SDIRK 2,\n\t"
+                  "7  - SDIRK 3,\n\t"
+                  "8  - EPIC (exponential euler)\n\t");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&alpha, "-a", "--alpha",
+                  "Alpha coefficient.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "Kappa coefficient offset.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (ode_solver_type < 1 || ode_solver_type > 9)
+   {
+      cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+      return 3;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
+   //    command-line parameter.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define the vector finite element space representing the current and the
+   //    initial temperature, u_ref.
+   H1_FECollection fe_coll(order, dim);
+   FiniteElementSpace fespace(mesh, &fe_coll);
+
+   int fe_size = fespace.GetTrueVSize();
+   cout << "Number of temperature unknowns: " << fe_size << endl;
+
+   GridFunction u_gf(&fespace);
+
+   // 5. Set the initial conditions for u. All boundaries are considered
+   //    natural.
+   FunctionCoefficient u_0(InitialTemperature);
+   u_gf.ProjectCoefficient(u_0);
+   Vector u;
+   u_gf.GetTrueDofs(u);
+
+   // 6. Initialize the conduction operator and the visualization.
+   ConductionOperator oper(fespace, alpha, kappa, u);
+
+   u_gf.SetFromTrueDofs(u);
+   {
+      ofstream omesh("ex16.mesh");
+      omesh.precision(precision);
+      mesh->Print(omesh);
+      ofstream osol("ex16-init.gf");
+      osol.precision(precision);
+      u_gf.Save(osol);
+   }
+
+   VisItDataCollection visit_dc("Example16", mesh);
+   visit_dc.RegisterField("temperature", &u_gf);
+   if (visit)
+   {
+      visit_dc.SetCycle(0);
+      visit_dc.SetTime(0.0);
+      visit_dc.Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << *mesh << u_gf;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 7. Define the ODE solver used for time integration.
+   double t = 0.0;
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // MFEM explicit methods
+      case 1: ode_solver = new ForwardEulerSolver; break;
+      case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 3: ode_solver = new RK3SSPSolver; break;
+      case 4: ode_solver = new RK4Solver; break;
+      // MFEM implicit L-stable methods
+      case 5: ode_solver = new BackwardEulerSolver; break;
+      case 6: ode_solver = new SDIRK23Solver(2); break;
+      case 7: ode_solver = new SDIRK33Solver; break;
+      // EPIC 
+      case 8: ode_solver = new EPI2();break;
+      case 9: ode_solver = new EPIRK4(); break;
+   }
+
+   // Initialize integrators
+   ode_solver->Init(oper);
+
+   // 8. Perform time-integration (looping over the time iterations, ti, with a
+   //    time-step dt).
+   cout << "Integrating the ODE ..." << endl;
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   /*ofstream out_jac_an("jacobian_an.txt");
+   ofstream out_jac_fd("jacobian_fd.txt");
+   oper.ExtractJacobians(u, out_jac_fd, out_jac_an);*/
+
+   bool last_step = false;
+   int ti;
+   for (ti = 1; !last_step; ti++)
+   {
+      double dt_real = min(dt, t_final - t);
+
+      // Note that since we are using the "one-step" mode of the SUNDIALS
+      // solvers, they will, generally, step over the final time and will not
+      // explicitly perform the interpolation to t_final as they do in the
+      // "normal" step mode.
+      ode_solver->Step(u, t, dt_real);
+
+      oper.UpdateStats();
+
+      last_step = (t >= t_final - 1e-8*dt);
+
+      if (last_step || (ti % vis_steps) == 0) {
+         cout << "step " << ti << ", t = " << t << endl;
+
+         u_gf.SetFromTrueDofs(u);
+         if (visualization) {
+            sout << "solution\n" << *mesh << u_gf << flush;
+         }
+
+         if (visit) {
+            visit_dc.SetCycle(ti);
+            visit_dc.SetTime(t);
+            visit_dc.Save();
+         }
+      }
+   }
+   tic_toc.Stop();
+   double comp_time = tic_toc.RealTime();
+   cout << "Done, " << comp_time << "s." << endl;
+
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -m ex16.mesh -g ex16-final.gf".
+   {
+      ofstream osol("ex16-final.gf");
+      osol.precision(precision);
+      u_gf.Save(osol);
+
+      ofstream ostats("ex16-stats.txt");
+      ostats << "time " << comp_time << endl;
+      oper.PrintStats(ostats);
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete mesh;
+
+   return 0;
+}
+
+ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al, double kap, const Vector &u)
+   : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL), dK(NULL), J_K(NULL), jac(NULL), z(height), u0(height),
+     nRhsMult(0), nSetJac(0), nJacMult(0), nImpSolve(0), nImpIter(0), nImpMult(0), nImpSet(0)
+{
+   const double rel_tol = 1e-8;
+
+   M = new BilinearForm(&fespace);
+   M->AddDomainIntegrator(new MassIntegrator());
+   M->Assemble();
+   M->FormSystemMatrix(ess_tdof_list, Mmat);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(rel_tol);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(50);
+   M_solver.SetPrintLevel(0);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(Mmat);
+
+   Jg_solver.SetRelTol(rel_tol);
+   Jg_solver.SetAbsTol(0.0);
+   Jg_solver.SetMaxIter(50);
+   Jg_solver.SetPrintLevel(0);
+   Jg_solver.SetPreconditioner(Jg_prec);
+
+   newton_solver.SetMaxIter(10);
+   newton_solver.SetRelTol(rel_tol);
+   newton_solver.SetPrintLevel(-1);
+   newton_solver.SetSolver(Jg_solver);
+   newton_solver.SetMaxIter(100);
+   newton_solver.iterative_mode = false;
+
+   alpha = al;
+   kappa = kap;
+}
+
+void ConductionOperator::UpdateStats()
+{
+   if (jac)
+   {
+      nJacMult += jac->GetnMult();
+   }
+}
+
+void ConductionOperator::PrintStats(ostream &out)
+{
+   out  << "nRhsMult " << nRhsMult << endl
+        << "nSetJac " << nSetJac << endl
+        << "nJacMult " << nJacMult  << endl
+        << "nImplicitSolve " << nImpSolve  << endl
+        << "nImplicitIter " << nImpIter << endl
+        << "nImplicitMult " << nImpMult << endl
+        << "nImplicitSet " << nImpSet << endl;
+}
+
+BilinearForm& ConductionOperator::GetKLambda(const Vector &u) const
+{
+   GridFunction conductivity_gf(&fespace);
+   conductivity_gf.SetFromTrueDofs(u);
+   for (int i = 0; i < conductivity_gf.Size(); i++)
+   {
+      conductivity_gf(i) = kappa + alpha*conductivity_gf(i);
+   }
+
+   GridFunctionCoefficient conductivity_coeff(&conductivity_gf);
+
+   delete K;
+   K = new BilinearForm(&fespace);
+   K->AddDomainIntegrator(new DiffusionIntegrator(conductivity_coeff));
+   K->Assemble();
+
+   return *K;
+}
+
+BilinearForm& ConductionOperator::GetdKLambda(const Vector &u) const
+{
+   GridFunction conductivity_gf(&fespace);
+   conductivity_gf.SetFromTrueDofs(u);
+   for (int i = 0; i < conductivity_gf.Size(); i++)
+   {
+      conductivity_gf(i) = kappa + alpha*conductivity_gf(i);
+   }
+
+   // Define diffusion form with conductivity = kappa(u0)
+   GridFunctionCoefficient conductivity_coeff(&conductivity_gf);
+
+   // Define advection form with velocity = grad kappa(u0)
+   GridFunction neg_cond_gf(conductivity_gf);
+   neg_cond_gf.Neg();
+   GradientGridFunctionCoefficient velocity_coeff(&neg_cond_gf);
+
+   delete dK;
+   dK = new BilinearForm(&fespace);
+
+   dK->AddDomainIntegrator(new DiffusionIntegrator(conductivity_coeff));
+   dK->AddDomainIntegrator(new MixedScalarWeakDivergenceIntegrator(velocity_coeff));
+   dK->Assemble();
+
+   return *dK;
+}
+
+void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
+{
+   // Compute:
+   //    du_dt = M^{-1}*-K(u)
+   // for du_dt
+   GetKLambda(u);
+   K->Mult(u, z);
+   z.Neg(); // z = -z
+   M_solver.Mult(z, du_dt);
+   nRhsMult++;
+}
+
+void ConductionOperator::ImplicitSolve(const double dt, const Vector &x, Vector &k)
+{
+   ImplicitSolveOperator imp_oper(this, &this->Mmat, dt, &x);
+   newton_solver.SetOperator(imp_oper);
+
+   Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
+   newton_solver.Mult(zero, k);
+   MFEM_VERIFY(newton_solver.GetConverged(), "Newton solver did not converge.");
+
+   nImpSolve++;
+   nImpMult += imp_oper.GetnMult();
+   nImpSet  += imp_oper.GetnSet();
+   nImpIter += newton_solver.GetNumIterations();
+}
+
+Operator &ConductionOperator::GetGradient(const Vector &u) const
+{
+   delete jac;
+   GetdKLambda(u);
+   jac = new JacobianOperator(dK, &M_solver);
+
+   nSetJac++;
+
+   return *jac;
+}
+
+ConductionOperator::~ConductionOperator()
+{
+   delete M;
+   delete K;
+   delete dK;
+   delete J_K;
+   delete jac;
+}
+
+ImplicitSolveOperator::ImplicitSolveOperator(ConductionOperator *oper_, const SparseMatrix* M_, double dt_, const Vector* x_):
+   Operator(oper_->Height()), oper(oper_), M(M_), dt(dt_), x(x_), u(height), z(height), Jg(NULL), nMult(0), nSet(0)
+{ }
+
+
+void ImplicitSolveOperator::Mult(const Vector& y, Vector& gy) const
+{
+   // Compute gy = g(y) = My + dt K(lambda(u)) u
+   // with u = x + dt y
+   add(*x, dt, y, u);
+   BilinearForm& K = oper->GetKLambda(u);
+   K.Mult(u, gy);
+
+   M->AddMult(y, gy);
+
+   nMult++;
+}
+
+Operator& ImplicitSolveOperator::GetGradient(const Vector &k) const
+{
+   add(*x, dt, k, u);
+
+   BilinearForm& dK = oper->GetdKLambda(u);
+   Array<int> ess_tdof_list;
+   SparseMatrix dK_mat;
+   dK.FormSystemMatrix(ess_tdof_list, dK_mat);
+
+   delete Jg;
+   Jg = Add(1.0, *M, dt, dK_mat);
+
+   nSet++;
+   return *Jg;
+}
+
+JacobianOperator::JacobianOperator(Operator* J_, Operator* M_solver_):
+   Operator(M_solver_->Height()), J(J_), M_solver(M_solver_), z(height), nMult(0)
+{ }
+
+void JacobianOperator::Mult(const Vector &v, Vector &Jv) const
+{
+   Vector temp(v);
+   J->Mult(v, z);
+   z.Neg(); // z = -z
+   M_solver->Mult(z, Jv);
+   nMult++;
+}
+
+
+void ConductionOperator::ExtractJacobians(const Vector& x, std::ostream &out, std::ostream &out2)
+{
+   int n = x.Size();
+
+   Vector e(n);
+   e = 0.0;
+
+   double eps = 1e-8;
+   Vector fx(n), fx_eps(n), x_eps(n);
+   Mult(x, fx);
+
+   DenseMatrix J(n);
+
+   for (int i = 0; i < n; i++)
+   {
+      e[i] = 1.0;
+      add(x, eps, e, x_eps);
+      Mult(x_eps, fx_eps);
+      fx_eps -= fx;
+      fx_eps /= eps;
+      J.SetCol(i, fx_eps);
+      e[i] = 0.0;
+   }
+
+   J.PrintMatlab(out);
+   GetGradient(x);
+   jac->ExtractJacobian(x, out2);
+}
+
+void JacobianOperator::ExtractJacobian(const Vector& x, std::ostream &out)
+{
+   int n = z.Size();
+
+   Vector e(n);
+   e= 0.0;
+
+   Vector J_i(n);
+   DenseMatrix J(n);
+
+   for (int i = 0; i < n; i++)
+   {
+      e[i] = 1.0;
+      Mult(e, J_i);
+      J.SetCol(i, J_i);
+      e[i] = 0.0;
+   }
+
+   J.PrintMatlab(out);
+}
+
+double InitialTemperature(const Vector &x)
+{
+   if (x.Norml2() < 0.5) { return 2.0; }
+   else                  { return 1.0; }
+}

examples/epic/ex16p.cpp

@@ -0,0 +1,494 @@
+//                         MFEM Example 16 - Parallel Version
+//                             SUNDIALS Modification
+//
+// Compile with: make ex16p
+//
+// Sample runs:
+//     mpirun -np 4 ex16p
+//     mpirun -np 4 ex16p -m ../../data/inline-tri.mesh
+//     mpirun -np 4 ex16p -m ../../data/disc-nurbs.mesh -tf 2
+//     mpirun -np 4 ex16p -s 12 -a 0.0 -k 1.0
+//     mpirun -np 4 ex16p -s 8 -a 1.0 -k 0.0 -dt 4e-6 -tf 2e-2 -vs 50
+//     mpirun -np 8 ex16p -s 9 -a 0.5 -k 0.5 -o 4 -dt 8e-6 -tf 2e-2 -vs 50
+//     mpirun -np 4 ex16p -s 10 -dt 2.0e-4 -tf 4.0e-2
+//     mpirun -np 16 ex16p -m ../../data/fichera-q2.mesh
+//     mpirun -np 16 ex16p -m ../../data/escher-p2.mesh
+//     mpirun -np 8 ex16p -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
+//     mpirun -np 4 ex16p -m ../../data/amr-quad.mesh -o 4 -rs 0 -rp 0
+//     mpirun -np 4 ex16p -m ../../data/amr-hex.mesh -o 2 -rs 0 -rp 0
+//
+// Description:  This example solves a time dependent nonlinear heat equation
+//               problem of the form du/dt = C(u), with a non-linear diffusion
+//               operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
+//
+//               We recommend viewing examples 2, 9 and 10 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+/** After spatial discretization, the conduction model can be written as:
+ *
+ *     du/dt = M^{-1}(-Ku)
+ *
+ *  where u is the vector representing the temperature, M is the mass matrix,
+ *  and K is the diffusion operator with diffusivity depending on u:
+ *  (\kappa + \alpha u).
+ *
+ *  Class ConductionOperator represents the right-hand side of the above ODE.
+ */
+class ConductionOperator : public TimeDependentOperator
+{
+protected:
+   ParFiniteElementSpace &fespace;
+   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
+
+   ParBilinearForm *M;
+   ParBilinearForm *K;
+
+   HypreParMatrix Mmat;
+   HypreParMatrix Kmat;
+   HypreParMatrix *T; // T = M + dt K
+   double current_dt;
+
+   CGSolver M_solver;    // Krylov solver for inverting the mass matrix M
+   HypreSmoother M_prec; // Preconditioner for the mass matrix M
+
+   CGSolver T_solver;    // Implicit solver for T = M + dt K
+   HypreSmoother T_prec; // Preconditioner for the implicit solver
+
+   double alpha, kappa;
+
+   mutable Vector z; // auxiliary vector
+
+public:
+   ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
+                      const Vector &u);
+
+   virtual void Mult(const Vector &u, Vector &du_dt) const;
+
+   /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
+       This is the only requirement for high-order SDIRK implicit integration.*/
+   virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
+
+   /** Setup the system (M + dt K) x = M b. This method is used by the implicit
+       SUNDIALS solvers. */
+   virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
+                                int jok, int *jcur, double gamma);
+
+   /** Solve the system (M + dt K) x = M b. This method is used by the implicit
+       SUNDIALS solvers. */
+   virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
+
+   /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
+   void SetParameters(const Vector &u);
+
+   virtual ~ConductionOperator();
+};
+
+double InitialTemperature(const Vector &x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int order = 2;
+   int ode_solver_type = 8; // Exponential Euler 
+   double t_final = 0.5;
+   double dt = 1.0e-2;
+   double alpha = 1.0e-2;
+   double kappa = 0.5;
+   bool visualization = true;
+   bool visit = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver:\n\t"
+                  "1  - Forward Euler,\n\t"
+                  "2  - RK2,\n\t"
+                  "3  - RK3 SSP,\n\t"
+                  "4  - RK4,\n\t"
+                  "5  - Backward Euler,\n\t"
+                  "6  - SDIRK 2,\n\t"
+                  "7  - SDIRK 3,\n\t"
+                  "8  - Exponential Euler,\n\t");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&alpha, "-a", "--alpha",
+                  "Alpha coefficient.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "Kappa coefficient offset.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      MPI_Finalize();
+      return 1;
+   }
+
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // check for vaild ODE solver option
+   if (ode_solver_type < 1 || ode_solver_type > 8)
+   {
+      if (myid == 0)
+      {
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+      }
+      MPI_Finalize();
+      return 1;
+   }
+
+   // 3. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle triangular, quadrilateral, tetrahedral and hexahedral meshes
+   //    with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define the vector finite element space representing the current and the
+   //    initial temperature, u_ref.
+   H1_FECollection fe_coll(order, dim);
+   ParFiniteElementSpace fespace(pmesh, &fe_coll);
+
+   int fe_size = fespace.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of temperature unknowns: " << fe_size << endl;
+   }
+
+   ParGridFunction u_gf(&fespace);
+
+   // 7. Set the initial conditions for u. All boundaries are considered
+   //    natural.
+   FunctionCoefficient u_0(InitialTemperature);
+   u_gf.ProjectCoefficient(u_0);
+   Vector u;
+   u_gf.GetTrueDofs(u);
+
+   // 8. Initialize the conduction operator and the VisIt visualization.
+   ConductionOperator oper(fespace, alpha, kappa, u);
+
+   u_gf.SetFromTrueDofs(u);
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "ex16-mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex16-init." << setfill('0') << setw(6) << myid;
+      ofstream omesh(mesh_name.str().c_str());
+      omesh.precision(precision);
+      pmesh->Print(omesh);
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u_gf.Save(osol);
+   }
+
+   VisItDataCollection visit_dc("Example16-Parallel", pmesh);
+   visit_dc.RegisterField("temperature", &u_gf);
+   if (visit)
+   {
+      visit_dc.SetCycle(0);
+      visit_dc.SetTime(0.0);
+      visit_dc.Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      sout << "parallel " << num_procs << " " << myid << endl;
+      int good = sout.good(), all_good;
+      MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm());
+      if (!all_good)
+      {
+         sout.close();
+         visualization = false;
+         if (myid == 0)
+         {
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+            cout << "GLVis visualization disabled.\n";
+         }
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << *pmesh << u_gf;
+         sout << "pause\n";
+         sout << flush;
+         if (myid == 0)
+         {
+            cout << "GLVis visualization paused."
+                 << " Press space (in the GLVis window) to resume it.\n";
+         }
+      }
+   }
+
+   // 9. Define the ODE solver used for time integration.
+   double t = 0.0;
+   ODESolver *ode_solver = NULL;
+   EPICSolver *epic_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // MFEM explicit methods
+      case 1: ode_solver = new ForwardEulerSolver; break;
+      case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 3: ode_solver = new RK3SSPSolver; break;
+      case 4: ode_solver = new RK4Solver; break;
+      // MFEM implicit L-stable methods
+      case 5: ode_solver = new BackwardEulerSolver; break;
+      case 6: ode_solver = new SDIRK23Solver(2); break;
+      case 7: ode_solver = new SDIRK33Solver; break;
+      // EPIC
+      case 8:
+         epic_solver = new EPICSolver();
+	 epic_solver->Init(oper);
+	 ode_solver = epic_solver;
+	 break;
+   }
+
+   // Initialize MFEM integrators
+   ode_solver->Init(oper);
+
+   // 10. Perform time-integration (looping over the time iterations, ti, with a
+   //     time-step dt).
+   if (myid == 0)
+   {
+      cout << "Integrating the ODE ..." << endl;
+   }
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   bool last_step = false;
+   for (int ti = 1; !last_step; ti++)
+   {
+      double dt_real = min(dt, t_final - t);
+
+      // Note that since we are using the "one-step" mode of the SUNDIALS
+      // solvers, they will, generally, step over the final time and will not
+      // explicitly perform the interpolation to t_final as they do in the
+      // "normal" step mode.
+
+      ode_solver->Step(u, t, dt_real);
+
+      last_step = (t >= t_final - 1e-8*dt);
+
+      if (last_step || (ti % vis_steps) == 0)
+      {
+         if (myid == 0)
+         {
+            cout << "step " << ti << ", t = " << t << endl;
+         }
+
+         u_gf.SetFromTrueDofs(u);
+         if (visualization)
+         {
+            sout << "parallel " << num_procs << " " << myid << "\n";
+            sout << "solution\n" << *pmesh << u_gf << flush;
+         }
+
+         if (visit)
+         {
+            visit_dc.SetCycle(ti);
+            visit_dc.SetTime(t);
+            visit_dc.Save();
+         }
+      }
+      oper.SetParameters(u);
+   }
+   tic_toc.Stop();
+   if (myid == 0)
+   {
+      cout << "Done, " << tic_toc.RealTime() << "s." << endl;
+   }
+
+   // 11. Save the final solution in parallel. This output can be viewed later
+   //     using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
+   {
+      ostringstream sol_name;
+      sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u_gf.Save(osol);
+   }
+
+   // 12. Free the used memory.
+   delete ode_solver;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
+                                       double kap, const Vector &u)
+   : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
+     T(NULL),
+     M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
+{
+   const double rel_tol = 1e-8;
+
+   M = new ParBilinearForm(&fespace);
+   M->AddDomainIntegrator(new MassIntegrator());
+   M->Assemble(0); // keep sparsity pattern of M and K the same
+   M->FormSystemMatrix(ess_tdof_list, Mmat);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(rel_tol);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(Mmat);
+
+   alpha = al;
+   kappa = kap;
+
+   T_solver.iterative_mode = false;
+   T_solver.SetRelTol(rel_tol);
+   T_solver.SetAbsTol(0.0);
+   T_solver.SetMaxIter(100);
+   T_solver.SetPrintLevel(0);
+   T_solver.SetPreconditioner(T_prec);
+
+   SetParameters(u);
+}
+
+void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
+{
+   // Compute:
+   //    du_dt = M^{-1}*-K(u)
+   // for du_dt
+   Kmat.Mult(u, z);
+   z.Neg(); // z = -z
+   M_solver.Mult(z, du_dt);
+}
+
+void ConductionOperator::ImplicitSolve(const double dt,
+                                       const Vector &u, Vector &du_dt)
+{
+   // Solve the equation:
+   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
+   // for du_dt
+   if (T) { delete T; }
+   T = Add(1.0, Mmat, dt, Kmat);
+   T_solver.SetOperator(*T);
+   Kmat.Mult(u, z);
+   z.Neg();
+   T_solver.Mult(z, du_dt);
+}
+
+int ConductionOperator::SUNImplicitSetup(const Vector &x,
+                                         const Vector &fx, int jok, int *jcur,
+                                         double gamma)
+{
+   // Setup the ODE Jacobian T = M + gamma K.
+   if (T) { delete T; }
+   T = Add(1.0, Mmat, gamma, Kmat);
+   T_solver.SetOperator(*T);
+   *jcur = 1;
+   return (0);
+}
+
+int ConductionOperator::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
+{
+   // Solve the system A x = z => (M - gamma K) x = M b.
+   Mmat.Mult(b, z);
+   T_solver.Mult(z, x);
+   return (0);
+}
+
+void ConductionOperator::SetParameters(const Vector &u)
+{
+   ParGridFunction u_alpha_gf(&fespace);
+   u_alpha_gf.SetFromTrueDofs(u);
+   for (int i = 0; i < u_alpha_gf.Size(); i++)
+   {
+      u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
+   }
+
+   delete K;
+   K = new ParBilinearForm(&fespace);
+
+   GridFunctionCoefficient u_coeff(&u_alpha_gf);
+
+   K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
+   K->Assemble(0); // keep sparsity pattern of M and K the same
+   K->FormSystemMatrix(ess_tdof_list, Kmat);
+}
+
+ConductionOperator::~ConductionOperator()
+{
+   delete T;
+   delete M;
+   delete K;
+}
+
+double InitialTemperature(const Vector &x)
+{
+   if (x.Norml2() < 0.5)
+   {
+      return 2.0;
+   }
+   else
+   {
+      return 1.0;
+   }
+}

examples/epic/makefile

@@ -0,0 +1,76 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/epic/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex16
+PAR_EXAMPLES = ex16p
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+ifeq ($(MFEM_USE_EPIC),NO)
+$(EXAMPLES):
+	$(error MFEM is not configured with EPIC)
+endif
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+SERIAL_NAME := Serial EPIC example
+PARALLEL_NAME := Parallel EPIC example
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME))
+%-test-seq: %
+	@$(call mfem-test,$<,, $(SERIAL_NAME))
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -f deformed.* velocity.* elastic_energy.*
+	@rm -f ex16.mesh ex16-mesh.* ex16-init.* ex16-final.* Example16*

examples/makefile

@@ -45,6 +45,9 @@ endif
 ifeq ($(MFEM_USE_HIOP),YES)
    SUBDIRS += hiop
 endif
+ifeq ($(MFEM_USE_EPIC),YES)
+   SUBDIRS += epic
+endif
 ifeq ($(MFEM_USE_PETSC),YES)
    SUBDIRS += petsc
 endif

general/version.cpp

@@ -82,6 +82,9 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_CUDA
       "MFEM_USE_CUDA\n"
 #endif
+#ifdef MFEM_USE_EPIC
+      "MFEM_USE_EPIC\n"
+#endif
 #ifdef MFEM_USE_EXCEPTIONS
       "MFEM_USE_EXCEPTIONS\n"
 #endif

linalg/CMakeLists.txt

@@ -83,6 +83,11 @@ if (MFEM_USE_SUNDIALS)
   list(APPEND HDRS sundials.hpp)
 endif()
 
+if (MFEM_USE_EPIC)
+  list(APPEND SRCS epic.cpp)
+  list(APPEND HDRS epic.hpp)
+endif()
+
 if (MFEM_USE_SUPERLU)
   list(APPEND SRCS superlu.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the

linalg/epic.cpp

@@ -0,0 +1,171 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "epic.hpp"
+
+#ifdef MFEM_USE_EPIC
+
+namespace mfem
+{
+
+EPICSolver::EPICSolver(bool exactJacobian_, EPICNumJacDelta delta)
+{
+   // Allocate an empty serial N_Vector
+   temp = N_VNewEmpty_Serial(0);
+   m[0] = 10;
+   m[1] = 10;
+   MFEM_VERIFY(temp, "error in N_VNewEmpty_Serial()");
+   exactJacobian = exactJacobian_;
+   Jtv = NULL;
+   Delta = delta;
+}
+
+#ifdef MFEM_USE_MPI
+EPICSolver::EPICSolver(MPI_Comm comm)
+{
+   m[0] = 10;
+   m[1] = 10;
+
+   // Allocate an empty vector
+   if (comm == MPI_COMM_NULL)
+   {
+      // Allocate an empty serial N_Vector
+      temp = N_VNewEmpty_Serial(0);
+      MFEM_VERIFY(temp, "error in N_VNewEmpty_Serial()");
+   }
+   else
+   {
+      // Allocate an empty parallel N_Vector
+      temp = N_VNewEmpty_Parallel(comm, 0, 0);  // calls MPI_Allreduce()
+      MFEM_VERIFY(temp, "error in N_VNewEmpty_Parallel()");
+   }
+}
+#endif
+
+int EPICSolver::RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
+{
+   // Get data from N_Vectors
+   const Vector mfem_y(y);
+   Vector mfem_ydot(ydot);
+   EPICSolver *self = static_cast<EPICSolver*>(user_data);
+
+   // Compute y' = f(t, y)
+   self->f->SetTime(t);
+   self->f->Mult(mfem_y, mfem_ydot);
+
+   // Return success
+   return 0;
+}
+
+
+int EPICSolver::Jacobian(N_Vector v, N_Vector Jv, realtype t, N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
+{
+    // Get data from N_Vectors
+   const Vector mfem_v(v);
+   Vector mfem_Jv(Jv);
+   EPICSolver *self = static_cast<EPICSolver*>(user_data);
+
+   // Compute J(t, y) v
+   self->Jtv->Mult(mfem_v, mfem_Jv);
+
+   return 0;
+}
+
+void EPICSolver::Init(TimeDependentOperator &f)
+{
+    ODESolver::Init(f);
+	    
+    long local_size = f.Height();
+    long global_size = 0;
+#ifdef MFEM_USE_MPI
+    if (Parallel())
+    {
+        MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
+                      NV_COMM_P(temp));
+    }
+#endif
+
+    Vector mfem_temp(local_size);
+    mfem_temp.ToNVector(temp, global_size);
+}
+
+EPI2::EPI2(bool exactJacobian, EPICNumJacDelta delta) : EPICSolver(exactJacobian, delta) {}
+
+void EPI2::Init(TimeDependentOperator &f)
+{
+    EPICSolver::Init(f);
+    long local_size = f.Height();
+    if (exactJacobian) {
+       integrator = new Epi2_KIOPS(EPICSolver::RHS, EPICSolver::Jacobian, this, 100, temp ,local_size);
+    } else {
+       integrator = new Epi2_KIOPS(EPICSolver::RHS, Delta, this, 100, temp ,local_size);
+    }
+
+}
+
+EPIRK4::EPIRK4(bool exactJacobian, EPICNumJacDelta delta) : EPICSolver(exactJacobian, delta) {}
+
+void EPIRK4::Init(TimeDependentOperator &f)
+{
+    EPICSolver::Init(f);
+    long local_size = f.Height();
+    if (exactJacobian) {
+       integrator = new EpiRK4SC_KIOPS(EPICSolver::RHS, EPICSolver::Jacobian, this, 100, temp ,local_size);
+    } else {
+       integrator = new EpiRK4SC_KIOPS(EPICSolver::RHS, Delta, this, 100, temp ,local_size);
+    }
+}
+
+void EPICSolver::Step(Vector &x, double &t, double &dt)
+{
+   if (!Parallel())
+   {
+      NV_DATA_S(temp) = x.GetData();
+      MFEM_VERIFY(NV_LENGTH_S(temp) == x.Size(), "");
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      NV_DATA_P(temp) = x.GetData();
+      MFEM_VERIFY(NV_LOCLENGTH_P(temp) == x.Size(), "");
+#endif
+   }
+
+   Jtv = &(this->f->GetGradient(x));
+}
+
+void EPI2::Step(Vector &x, double &t, double &dt)
+{
+    EPICSolver::Step(x, t, dt);
+    integrator->Integrate(dt, t, t+dt, 0, temp, 1e-10, m);
+    t += dt;
+}
+
+void EPIRK4::Step(Vector &x, double &t, double &dt)
+{
+    EPICSolver::Step(x, t, dt);
+    integrator->Integrate(dt, t, t+dt, 0, temp, 1e-10, m);
+    t += dt;
+}
+
+EPI2::~EPI2()
+{
+    delete integrator;
+}
+
+EPIRK4::~EPIRK4()
+{
+    delete integrator;
+}
+
+}
+
+#endif

linalg/epic.hpp

@@ -0,0 +1,98 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_EPIC
+#define MFEM_EPIC
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_EPIC
+
+// SUNDIALS vectors
+#include <nvector/nvector_serial.h>
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#include <nvector/nvector_parallel.h>
+#endif
+
+#include "ode.hpp"
+#include "solvers.hpp"
+#include <Epic.h>
+
+namespace mfem
+{
+
+typedef void (*JacobianFun)(const realtype t, const Vector &y, const Vector& v, Vector& Jv, void* user_data);
+// ---------------------------------------------------------------------------
+// Interface to the EPIC library -- exponential methods
+// ---------------------------------------------------------------------------
+
+class EPICSolver : public ODESolver
+{
+protected:
+    EPICNumJacDelta Delta;
+    Operator* Jtv;
+    N_Vector temp;
+    int m[2];
+
+    bool exactJacobian;
+
+    #ifdef MFEM_USE_MPI
+    bool Parallel() const
+    {
+        return (N_VGetVectorID(temp) != SUNDIALS_NVEC_SERIAL);
+    }
+    #else
+    bool Parallel() const { return false; }
+    #endif
+
+public:
+    EPICSolver(bool exactJacobian, EPICNumJacDelta delta=&DefaultDelta);
+    EPICSolver(MPI_Comm comm);
+
+    static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
+    static int Jacobian(N_Vector v, N_Vector Jv, realtype t,
+                         N_Vector y, N_Vector fy, void *user_data, N_Vector tmp);
+    virtual void Init(TimeDependentOperator &f);
+    virtual void Step(Vector &x, double &t, double &dt);
+
+    virtual ~EPICSolver() {}
+};
+
+class EPI2 : public EPICSolver
+{
+protected:
+    Epi2_KIOPS* integrator;
+public:
+    EPI2(bool exactJacobian=true, EPICNumJacDelta delta=&DefaultDelta);
+    virtual void Init(TimeDependentOperator &f);
+    virtual void Step(Vector &x, double &t, double &dt);
+
+    virtual ~EPI2();
+};
+
+class EPIRK4 : public EPICSolver
+{
+protected:
+    EpiRK4SC_KIOPS* integrator;
+public:
+    EPIRK4(bool exactJacobian=true, EPICNumJacDelta delta=&DefaultDelta);
+    virtual void Init(TimeDependentOperator &f);
+    virtual void Step(Vector &x, double &t, double &dt);
+
+    virtual ~EPIRK4();
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_EPIC
+
+#endif // MFEM_EPIC

linalg/linalg.hpp

@@ -40,6 +40,10 @@
 #include "sundials.hpp"
 #endif
 
+#ifdef MFEM_USE_EPIC
+#include "epic.hpp"
+#endif
+
 #ifdef MFEM_USE_HIOP
 #include "hiop.hpp"
 #endif

linalg/sundials.cpp

@@ -425,6 +425,26 @@ static SUNMatrix_ID MatGetID(SUNMatrix)
    return (SUNMATRIX_CUSTOM);
 }
 
+//BEGIN WORKAROUND CODE
+static SUNMatrix MatClone(SUNMatrix A)
+{
+   SUNMatrix B = SUNMatNewEmpty();
+   MFEM_VERIFY(B, "error in MatClone()");
+
+   B->content = A->content;
+
+   int flag = SUNMatCopyOps(A, B);
+   MFEM_VERIFY(flag == SUNMAT_SUCCESS, "error in MatClone()");
+
+   return B;
+}
+
+static int MatCopy(SUNMatrix, SUNMatrix)
+{
+   return 0;
+}
+// END WORKAROUN CODE
+
 static void MatDestroy(SUNMatrix A)
 {
    if (A->content) { A->content = NULL; }
@@ -1493,7 +1513,12 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
    MFEM_VERIFY(M, "error in SUNMatNewEmpty()");
 
    M->content      = this;
-   M->ops->getid   = SUNMatGetID;
+// BEGIN WORKAROUND CODE
+   // M->ops->getid   = SUNMatGetID;
+   M->ops->getid = MatGetID;
+   M->ops->clone = MatClone;
+   M->ops->copy = MatCopy;
+// END WORKAROUND CODE
    M->ops->matvec  = ARKStepSolver::MassMult1;
    M->ops->destroy = MatDestroy;
 
@@ -1537,6 +1562,12 @@ void ARKStepSolver::SetSStolerances(double reltol, double abstol)
    MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSStolerances()");
 }
 
+void ARKStepSolver::SetMaxIter(int iterations)
+{
+  flag = ARKStepSetMaxNonlinIters(sundials_mem, iterations);
+  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMaxNonlinIters()");
+}
+
 void ARKStepSolver::SetMaxStep(double dt_max)
 {
    flag = ARKStepSetMaxStep(sundials_mem, dt_max);

linalg/sundials.hpp

@@ -666,6 +666,9 @@ public:
    /// Set the scalar relative and scalar absolute tolerances.
    void SetSStolerances(double reltol, double abstol);
 
+   /// Set the maximum number of nonlinear iterations per RK stage
+   void SetMaxIter(int iterations);
+
    /// Set the maximum time step.
    void SetMaxStep(double dt_max);
 

makefile

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -119,7 +119,7 @@ $(if $(word 2,$(SRC)),$(error Spaces in SRC = "$(SRC)" are not supported))
 MFEM_GIT_STRING = $(shell [ -d $(MFEM_DIR)/.git ] && git -C $(MFEM_DIR) \
    describe --all --long --abbrev=40 --dirty --always 2> /dev/null)
 
-EXAMPLE_SUBDIRS = amgx ginkgo hiop petsc pumi sundials superlu
+EXAMPLE_SUBDIRS = amgx ginkgo hiop petsc pumi sundials epic superlu
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
@@ -272,7 +272,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS EPIC MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
  OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
 
@@ -333,7 +333,7 @@ endif
 MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_METIS MFEM_USE_METIS_5 MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB\
  MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_OPENMP\
- MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
+ MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS MFEM_USE_EPIC\
  MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
@@ -650,6 +650,7 @@ status info:
 	$(info MFEM_USE_MEMALLOC      = $(MFEM_USE_MEMALLOC))
 	$(info MFEM_TIMER_TYPE        = $(MFEM_TIMER_TYPE))
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
+	$(info MFEM_USE_EPIC          = $(MFEM_USE_EPIC))
 	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))