compile error fixed

linalg/dual.hpp

@@ -403,8 +403,105 @@ template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
 dual<value_type, gradient_type> pow(dual<value_type, gradient_type> a, real_t b)
 {
    using std::pow;
-   value_type value = pow(a.value, b);
-   return {value, value * a.gradient * b / a.value};
+   return {pow(a.value, b), b*pow(a.value, b-1) * a.gradient };
+}
+
+/** @brief implementation of max of two dual numbers */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> max(dual<value_type, gradient_type> a,
+                                    dual<value_type, gradient_type> b)
+{
+   using std::max;
+   if (a.value > b.value)
+   {
+      return a;
+   }
+   else if (a.value < b.value)
+   {
+      return b;
+   }
+   else { return (a+b)*0.5; } // subgradient at the kink
+}
+
+/** @brief implementation of max of a dual number and a non-dual number */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> max(real_t a, dual<value_type, gradient_type> b)
+{
+   using std::max;
+   if (a > b.value)
+   {
+      return {a, {}};
+   }
+   else if (a < b.value)
+   {
+      return b;
+   }
+   else { return {a, b.gradient*0.5}; } // subgradient at the kink
+}
+
+/** @brief implementation of max of two non-dual numbers */
+template <typename value_type > MFEM_HOST_DEVICE
+value_type max(value_type a, value_type b)
+{
+   using std::pow;
+   return max(a, b);
+}
+
+/** @brief implementation of max of a dual number and a non-dual number */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> max(dual<value_type, gradient_type> a, real_t b)
+{
+   using std::max;
+   return max(b, a);
+}
+
+/** @brief implementation of min of two dual numbers */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> min(dual<value_type, gradient_type> a,
+                                    dual<value_type, gradient_type> b)
+{
+   using std::max;
+   if (a.value < b.value)
+   {
+      return a;
+   }
+   else if (a.value > b.value)
+   {
+      return b;
+   }
+   else { return (a+b)*0.5; } // subgradient at the kink
+}
+
+/** @brief implementation of min of a dual number and a non-dual number */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> min(real_t a, dual<value_type, gradient_type> b)
+{
+   using std::max;
+   if (a < b.value)
+   {
+      return {a, {}};
+   }
+   else if (a > b.value)
+   {
+      return b;
+   }
+   else { return {a, b.gradient*0.5}; } // subgradient at the kink
+}
+
+/** @brief implementation of min of two non-dual numbers */
+template <typename value_type > MFEM_HOST_DEVICE
+value_type min(value_type a, value_type b)
+{
+   using std::pow;
+   return min(a, b);
+}
+
+/** @brief implementation of min of a dual number and a non-dual number */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> min(dual<value_type, gradient_type> a, real_t b)
+{
+   using std::max;
+   return min(b, a);
 }
 
 /** @brief overload of operator<< for `dual` to work with work with standard output streams */

miniapps/autodiff/CMakeLists.txt

@@ -9,12 +9,22 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-list(APPEND SEQADIFF_COMMON_SOURCES)
+list(APPEND SEQADIFF_COMMON_SOURCES
+  ad_native.cpp
+  logger.cpp
+  pg.cpp
+)
 
 list(APPEND SEQADIFF_COMMON_HEADERS
 	tadvector.hpp
 	taddensemat.hpp
-	admfem.hpp)
+	admfem.hpp
+  ad_intg.hpp
+  ad_native.hpp
+  logger.hpp
+  pg.hpp
+  tools.hpp
+)
 
 convert_filenames_to_full_paths(SEQADIFF_COMMON_SOURCES)
 convert_filenames_to_full_paths(SEQADIFF_COMMON_HEADERS)
@@ -23,34 +33,15 @@ set(SEQADIFF_COMMON_FILES
 	EXTRA_SOURCES ${SEQADIFF_COMMON_SOURCES}
 	EXTRA_HEADERS ${SEQADIFF_COMMON_HEADERS})
 
-add_mfem_miniapp(seqadiff
-	MAIN seq_example.cpp
-	${SEQADIFF_COMMON_FILES}
-	LIBRARIES mfem)
-
-add_mfem_miniapp(seqtest
-	MAIN seq_test.cpp
-	${SEQADIFF_COMMON_FILES}
-        LIBRARIES mfem)
-
-if(MFEM_USE_MPI)
-
-list(APPEND PARADIFF_COMMON_SOURCES)
-list(APPEND PARADIFF_COMMON_HEADERS)
-
-convert_filenames_to_full_paths(PARADIFF_COMMON_SOURCES)
-convert_filenames_to_full_paths(PARADIFF_COMMON_HEADERS)
-
-set(PARADIFF_COMMON_FILES
-	EXTRA_SOURCES ${PARADIFF_COMMON_SOURCES} ${SEQADIFF_COMMON_SOURCES}
-	EXTRA_HEADERS ${PARADIFF_COMMON_HEADERS} ${SEQADIFF_COMMON_HEADERS})
-
-# message(STATUS "PARADIFF_COMMON_FILES: ${PARADIFF_COMMON_FILES}")
-# message(STATUS "SEQADIFF_COMMON_FILES: ${SEQADIFF_COMMON_FILES}")
-
-add_mfem_miniapp(paradiff
-     MAIN par_example.cpp
-     ${PARADIFF_COMMON_FILES}
-     LIBRARIES mfem)
-
+add_mfem_miniapp(ad_ex0 MAIN ad_ex0.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
+add_mfem_miniapp(ad_ex1 MAIN ad_ex1.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
+add_mfem_miniapp(ad_ex2 MAIN ad_ex2.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
+add_mfem_miniapp(ad_ex3 MAIN ad_ex3.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
+if (MFEM_USE_MUMPS)
+  add_mfem_miniapp(ad_ex4 MAIN ad_ex4.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
+  add_mfem_miniapp(ad_ex5 MAIN ad_ex5.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
 endif ()
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(ad_ex6 MAIN ad_ex6.cpp ${SEQADIFF_COMMON_FILES} LIBRARIES mfem)
+endif ()
+

miniapps/autodiff/_ad_intg.hpp

@@ -0,0 +1,697 @@
+// IWYU pragma: private, include "ad_intg.hpp"
+// -----------------------------------------
+/// Templated AD (block) nonlinear form integrators implementations
+#pragma once
+#include "ad_intg.hpp"
+
+namespace mfem
+{
+
+template <ADEval mode>
+inline int ADNonlinearFormIntegrator<mode>::InitInputShapes(
+   const FiniteElement &el,
+   ElementTransformation &Tr,
+   DenseMatrix &shapes)
+{
+   const int sdim = Tr.GetSpaceDim();
+   const int dim = el.GetDim();
+   int idx[static_cast<int>(ADEval::NUMOPT)];
+   idx[0] = 0;
+   idx[1] = idx[0] + (hasFlag(mode, ADEval::QVALUE) ? 1 : 0);
+   idx[2] = idx[1] +  (hasFlag(mode, ADEval::VALUE)
+                       ? hasFlag(mode, ADEval::VECFE)
+                       ? dim // if vector-FE
+                       : 1 // if scalar-FE
+                       : 0); // no value
+   idx[3] = idx[2] + (hasFlag(mode, ADEval::GRAD) ? sdim : 0);
+   idx[4] = idx[3] + (hasFlag(mode, ADEval::DIV) ? 1 : 0);
+   idx[5] = idx[4] + (hasFlag(mode, ADEval::CURL) ? el.GetCurlDim() : 0);
+   const int shapedim = idx[5];
+   const int dof = el.GetDof();
+   shapes.SetSize(dof, shapedim);
+
+   if constexpr (hasFlag(mode, ADEval::QVALUE)) { shapes.SetCol(idx[0], 0.0); }
+
+   if constexpr (hasFlag(mode, ADEval::VALUE))
+   {
+      if constexpr (hasFlag(mode, ADEval::VECFE)) { vshape.UseExternalData(shapes.GetData() + dof*idx[1], dof, dim); }
+      else { shapes.GetColumnReference(idx[1], shape); }
+   }
+
+   if constexpr (hasFlag(mode, ADEval::GRAD))
+   {
+      gshape.UseExternalData(shapes.GetData() + dof*idx[2],
+                             dof, sdim);
+   }
+   if constexpr (hasFlag(mode, ADEval::DIV))
+   {
+      shapes.GetColumnReference(idx[3], divshape);
+   }
+
+   if constexpr (hasFlag(mode, ADEval::CURL))
+   {
+      curlshape.UseExternalData(shapes.GetData() + dof*idx[4],
+                                dof, el.GetCurlDim());
+   }
+
+   return shapedim;
+}
+
+template <ADEval mode>
+inline void ADNonlinearFormIntegrator<mode>::CalcInputShapes(
+   const FiniteElement &el,
+   ElementTransformation &Tr,
+   const IntegrationPoint &ip,
+   DenseMatrix &allshapes)
+{
+   // Get quadrature value
+   // ip should be from the same integration rule with base quadrature
+   if constexpr (hasFlag(mode, ADEval::QVALUE)) { allshapes.SetCol(0, 0.0); allshapes(ip.index, 0) = 1.0; }
+
+   // Get value shape
+   if constexpr (hasFlag(mode, ADEval::VALUE))
+   {
+      if constexpr (hasFlag(mode, ADEval::VECFE)) { el.CalcPhysVShape(Tr, vshape); }
+      else { el.CalcPhysShape(Tr, shape); }
+   }
+
+   // Get gradient shape
+   if constexpr (hasFlag(mode, ADEval::GRAD)) { el.CalcPhysDShape(Tr, gshape); }
+
+   // Get divergence shape
+   if constexpr (hasFlag(mode, ADEval::DIV))
+   {
+      if constexpr (hasFlag(mode, ADEval::GRAD))
+      {
+         gshape.GetRowSums(divshape);
+      }
+      else
+      {
+         el.CalcPhysDivShape(Tr, divshape);
+      }
+   }
+
+   // Get divergence shape
+   if constexpr (hasFlag(mode, ADEval::CURL)) { el.CalcPhysCurlShape(Tr, curlshape); }
+}
+
+/// Perform the local action of the NonlinearFormIntegrator
+template <ADEval mode>
+real_t ADNonlinearFormIntegrator<mode>::GetElementEnergy(
+   const FiniteElement &el,
+   ElementTransformation &Tr,
+   const Vector &elfun)
+{
+   const int dof = el.GetDof();
+   const int vdim = elfun.Size() / dof;
+   MFEM_ASSERT(vdim == 1 ? true : hasFlag(mode, ADEval::VECTOR),
+               "ADNonlinearFormIntegrator: "
+               "vdim must be 1 or the mode must be VECTOR");
+
+   real_t energy = 0.0;
+
+   int shapedim = InitInputShapes(el, Tr, allshapes);
+   x.SetSize(f.n_input);
+   if constexpr (hasFlag(mode, ADEval::VECTOR))
+   {
+      elfun_matview.UseExternalData(const_cast<real_t*>(elfun.GetData()),
+                                    dof, vdim);
+      xmat.UseExternalData(x.GetData(), shapedim, vdim);
+   }
+
+   const IntegrationRule * ir = GetIntegrationRule(el, Tr);
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+
+      CalcInputShapes(el, Tr, ip, allshapes);
+
+      if constexpr (hasFlag(mode, ADEval::VECTOR))
+      {
+         MultAtB(allshapes, elfun_matview, xmat);
+      }
+      else
+      {
+         allshapes.MultTranspose(elfun, x);
+      }
+      energy += f(x, Tr, ip)*Tr.Weight()*ip.weight;
+   }
+   return energy;
+}
+
+/// Compute the local <grad f, v>
+template <ADEval mode>
+void ADNonlinearFormIntegrator<mode>::AssembleElementVector(
+   const FiniteElement &el,
+   ElementTransformation &Tr,
+   const Vector &elfun, Vector &elvect)
+{
+   const int dof = el.GetDof();
+   const int vdim = elfun.Size() / dof;
+   MFEM_ASSERT(vdim == 1 ? true : hasFlag(mode, ADEval::VECTOR),
+               "ADNonlinearFormIntegrator: "
+               "vdim must be 1 or the mode must be VECTOR");
+
+   real_t w;
+   elvect.SetSize(dof*vdim);
+   elvect = 0.0;
+
+   x.SetSize(f.n_input);
+   jac.SetSize(f.n_input);
+
+   int shapedim = InitInputShapes(el, Tr, allshapes);
+
+   if constexpr (hasFlag(mode, ADEval::VECTOR))
+   {
+      elfun_matview.UseExternalData(const_cast<real_t*>(elfun.GetData()),
+                                    dof, vdim);
+      elvectmat.UseExternalData(elvect.GetData(), dof, vdim);
+      xmat.UseExternalData(x.GetData(), shapedim, vdim);
+      jacMat.UseExternalData(jac.GetData(), shapedim, vdim);
+   }
+
+   const IntegrationRule * ir = GetIntegrationRule(el, Tr);
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      w = ip.weight * Tr.Weight();
+
+      CalcInputShapes(el, Tr, ip, allshapes);
+
+      // Convert dof to x = [[value, grad], [value, grad], ...]
+      if constexpr (hasFlag(mode, ADEval::VECTOR)) { MultAtB(allshapes, elfun_matview, xmat); }
+      else { allshapes.MultTranspose(elfun, x); }
+
+      f.Gradient(x, Tr, ip, jac);
+      jac *= w;
+
+      if constexpr (hasFlag(mode, ADEval::VECTOR))
+      {
+         AddMult(allshapes, jacMat, elvectmat);
+      }
+      else
+      {
+         allshapes.AddMult(jac, elvect);
+      }
+   }
+}
+
+/// Assemble the local <H_f(x)(u), v>
+template <ADEval mode>
+void ADNonlinearFormIntegrator<mode>::AssembleElementGrad(
+   const FiniteElement &el,
+   ElementTransformation &Tr,
+   const Vector &elfun, DenseMatrix &elmat)
+{
+   const int dof = el.GetDof();
+   const int vdim = elfun.Size() / dof;
+   MFEM_ASSERT(vdim == 1 ? true : hasFlag(mode, ADEval::VECTOR),
+               "ADNonlinearFormIntegrator: "
+               "vdim must be 1 or the mode must be VECTOR");
+
+   real_t w;
+   elmat.SetSize(dof*vdim);
+   elmat = 0.0;
+
+   int shapedim = InitInputShapes(el, Tr, allshapes);
+   MFEM_ASSERT(shapedim*vdim == f.n_input,
+               "ADNonlinearFormIntegrator: "
+               "shapedim*vdim must match n_input");
+
+   x.SetSize(f.n_input);
+   H.SetSize(f.n_input);
+   Hx.SetSize(dof, shapedim*vdim*vdim);
+
+
+   if constexpr (hasFlag(mode, ADEval::VECTOR))
+   {
+      elfun_matview.UseExternalData(const_cast<real_t*>(elfun.GetData()),
+                                    dof, vdim);
+      xmat.UseExternalData(x.GetData(), shapedim, vdim);
+      partelmat.SetSize(dof, dof);
+      Hs.UseExternalData(H.GetData(), shapedim, vdim*shapedim*vdim);
+   }
+
+   const IntegrationRule * ir = GetIntegrationRule(el, Tr);
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      w = ip.weight * Tr.Weight();
+      CalcInputShapes(el, Tr, ip, allshapes);
+
+      // Convert dof to x = [[value, grad], [value, grad], ...]
+      if constexpr (hasFlag(mode, ADEval::VECTOR)) { MultAtB(allshapes, elfun_matview, xmat); }
+      else { allshapes.MultTranspose(elfun, x); }
+
+      f.Hessian(x, Tr, ip, H);
+      H *= w;
+
+      if constexpr (hasFlag(mode, ADEval::VECTOR))
+      {
+         Mult(allshapes, Hs, Hx);
+         const int nel = shapedim*dof;
+         for (int c=0; c<vdim; c++)
+         {
+            for (int r=0; r<=c; r++)
+            {
+               Hxsub.UseExternalData(Hx.GetData() + (c*vdim + r)*nel, dof, shapedim);
+               MultABt(allshapes, Hxsub, partelmat);
+               elmat.AddSubMatrix(c*dof, r*dof, partelmat);
+               if (c != r)
+               {
+                  elmat.AddSubMatrix(r*dof, c*dof, partelmat);
+               }
+            }
+         }
+      }
+      else
+      {
+         Mult(allshapes, H, Hx);
+         AddMultABt(allshapes, Hx, elmat);
+      }
+   }
+}
+
+/// @brief Perform the local action of the NonlinearFormIntegrator resulting
+/// from a face integral term.
+template <ADEval mode>
+void ADNonlinearFormIntegrator<mode>::AssembleFaceVector(
+   const FiniteElement &el1,
+   const FiniteElement &el2,
+   FaceElementTransformations &Tr,
+   const Vector &elfun, Vector &elvect)
+{
+   MFEM_ABORT("ADNonlinearFormIntegrator::AssembleFaceVector: "
+              "This method is not implemented.");
+}
+
+
+/// @brief Assemble the local action of the gradient of the
+/// NonlinearFormIntegrator resulting from a face integral term.
+template <ADEval mode>
+void ADNonlinearFormIntegrator<mode>::AssembleFaceGrad(
+   const FiniteElement &el1,
+   const FiniteElement &el2,
+   FaceElementTransformations &Tr,
+   const Vector &elfun, DenseMatrix &elmat)
+{
+   MFEM_ABORT("ADNonlinearFormIntegrator::AssembleFaceVector: "
+              "This method is not implemented.");
+}
+
+template <ADEval... modes>
+inline std::array<int, sizeof...(modes)>
+ADBlockNonlinearFormIntegrator<modes...>::InitInputShapes(
+   const Array<const FiniteElement *>& els,
+   ElementTransformation &Tr,
+   std::vector<DenseMatrix> &shapes)
+{
+   MFEM_ASSERT(els.Size() == numSpaces,
+               "ADBlockNonlinearFormIntegrator: "
+               "el.Size()=" << els.Size() << " must match numSpaces=" << numSpaces);
+   const int sdim = Tr.GetSpaceDim();
+   std::array<int, sizeof...(modes)> shapedims{};
+
+   _constexpr_for([&](auto i)
+   {
+      constexpr auto mode = modes_arr[i];
+      const FiniteElement &el = *els[i];
+      const int sdim = Tr.GetSpaceDim();
+      const int dim = el.GetDim();
+      int idx[static_cast<int>(ADEval::NUMOPT)];
+      idx[0] = 0;
+      idx[1] = idx[0] + (hasFlag(modes_arr[i], ADEval::QVALUE) ? 1 : 0);
+      idx[2] = idx[1] +  (hasFlag(modes_arr[i], ADEval::VALUE)
+                          ? hasFlag(modes_arr[i], ADEval::VECFE)
+                          ? dim // if vector-FE
+                          : 1 // if scalar-FE
+                          : 0); // no value
+      idx[3] = idx[2] + (hasFlag(modes_arr[i], ADEval::GRAD) ? sdim : 0);
+      idx[4] = idx[3] + (hasFlag(modes_arr[i], ADEval::DIV) ? 1 : 0);
+      idx[5] = idx[4] + (hasFlag(modes_arr[i], ADEval::CURL) ? el.GetCurlDim() : 0);
+      const int shapedim = idx[5];
+      const int dof = el.GetDof();
+      shapes[i].SetSize(dof, shapedim);
+
+      if constexpr (hasFlag(mode, ADEval::QVALUE)) { shapes[i].SetCol(idx[0], 0.0); }
+
+      if constexpr (hasFlag(mode, ADEval::VALUE))
+      {
+         if constexpr (hasFlag(mode, ADEval::VECFE)) { vshape[i].UseExternalData(shapes[i].GetData() + dof*idx[1], dof, dim); }
+         else { shapes[i].GetColumnReference(idx[1], shape[i]); }
+      }
+
+      if constexpr (hasFlag(mode, ADEval::GRAD))
+      {
+         gshape[i].UseExternalData(shapes[i].GetData() + dof*idx[2],
+                                   dof, sdim);
+      }
+      if constexpr (hasFlag(mode, ADEval::DIV))
+      {
+         shapes[i].GetColumnReference(idx[3], divshape[i]);
+      }
+
+      if constexpr (hasFlag(mode, ADEval::CURL))
+      {
+         curlshape[i].UseExternalData(shapes[i].GetData() + dof*idx[4],
+                                      dof, el.GetCurlDim());
+      }
+      shapedims[i] = shapedim;
+   }, std::make_index_sequence<sizeof...(modes)> {});
+   return shapedims;
+}
+template <ADEval... modes>
+inline void
+ADBlockNonlinearFormIntegrator<modes...>::CalcInputShapes(
+   const Array<const FiniteElement *>& els,
+   ElementTransformation &Tr,
+   const IntegrationPoint &ip,
+   std::vector<DenseMatrix> &allshapes)
+{
+   _constexpr_for([&](auto i)
+   {
+      const auto&el = *els[i];
+      constexpr auto mode = modes_arr[i];
+      // Get quadrature value
+      // ip should be from the same integration rule with base quadrature
+      if constexpr (hasFlag(mode, ADEval::QVALUE)) { allshapes[i].SetCol(0, 0.0); allshapes[i](ip.index, 0) = 1.0; }
+
+      // Get value shape
+      if constexpr (hasFlag(mode, ADEval::VALUE))
+      {
+         if constexpr (hasFlag(mode, ADEval::VECFE)) { el.CalcPhysVShape(Tr, vshape[i]); }
+         else { el.CalcPhysShape(Tr, shape[i]); }
+      }
+
+      // Get gradient shape
+      if constexpr (hasFlag(mode, ADEval::GRAD)) { el.CalcPhysDShape(Tr, gshape[i]); }
+
+      // Get divergence shape
+      if constexpr (hasFlag(mode, ADEval::DIV))
+      {
+         if constexpr (hasFlag(mode, ADEval::GRAD))
+         {
+            gshape[i].GetRowSums(divshape[i]);
+         }
+         else
+         {
+            el.CalcPhysDivShape(Tr, divshape[i]);
+         }
+      }
+
+      // Get divergence shape
+      if constexpr (hasFlag(mode, ADEval::CURL)) { el.CalcPhysCurlShape(Tr, curlshape[i]); }
+   }, std::make_index_sequence<sizeof...(modes)> {});
+}
+
+/// Compute the local energy
+template <ADEval... modes>
+real_t ADBlockNonlinearFormIntegrator<modes...>::GetElementEnergy(
+   const Array<const FiniteElement *> &el,
+   ElementTransformation &Tr,
+   const Array<const Vector*> &elfun)
+{
+   MFEM_ASSERT(el.Size() == numSpaces,
+               "ADBlockNonlinearFormIntegrator: "
+               "el.Size()=" << el.Size() << " must match numSpaces=" << numSpaces);
+   std::array<int, numSpaces> dof{};
+   std::array<int, numSpaces> order{};
+   for (int i=0; i<numSpaces; i++)
+   {
+      dof[i] = el[i]->GetDof();
+      order[i] = el[i]->GetOrder();
+      vdim[i] = elfun[i]->Size() / dof[i];
+
+      MFEM_ASSERT(vdim[i] == 1 ? true : hasFlag(modes_arr[i], ADEval::VECTOR),
+                  "ADNonlinearFormIntegrator: "
+                  "vdim must be 1 or the mode must be VECTOR");
+   }
+
+   real_t energy = 0.0;
+
+   std::array<int, numSpaces> shapedim(InitInputShapes(el, Tr, allshapes));
+   x.SetSize(f.n_input);
+   int x_idx = 0;
+   _constexpr_for([&](auto vi)
+   {
+      xvar[vi].MakeRef(x, x_idx, shapedim[vi]*vdim[vi]);
+      x_idx += shapedim[vi]*vdim[vi];
+      if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+      {
+         elfun_matview[vi].UseExternalData(const_cast<real_t*>(elfun[vi]->GetData()),
+                                           dof[vi], vdim[vi]);
+         xmat[vi].UseExternalData(xvar[vi].GetData(), shapedim[vi], vdim[vi]);
+      }
+   }, std::make_index_sequence<sizeof...(modes)> {});
+
+   const IntegrationRule * ir = GetIntegrationRule(el, Tr);
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+
+      CalcInputShapes(el, Tr, ip, allshapes);
+
+      _constexpr_for([&](auto vi)
+      {
+         if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+         {
+            MultAtB(allshapes[vi], elfun_matview[vi], xmat[vi]);
+         }
+         else
+         {
+            allshapes[vi].MultTranspose(*elfun[vi], xvar[vi]);
+         }
+      }, std::make_index_sequence<sizeof...(modes)> {});
+      energy += f(x, Tr, ip)*Tr.Weight()*ip.weight;
+   }
+   return energy;
+}
+
+/// Perform the local action of the NonlinearFormIntegrator
+template <ADEval... modes>
+void ADBlockNonlinearFormIntegrator<modes...>::AssembleElementVector(
+   const Array<const FiniteElement *>&el,
+   ElementTransformation &Tr,
+   const Array<const Vector *>&elfun,
+   const Array<Vector *>&elvect)
+{
+   MFEM_ASSERT(el.Size() == numSpaces,
+               "ADBlockNonlinearFormIntegrator: "
+               "el.Size()=" << el.Size() << " must match numSpaces=" << numSpaces);
+   std::array<int, numSpaces> dof{};
+   std::array<int, numSpaces> order{};
+   for (int i=0; i<numSpaces; i++)
+   {
+      dof[i] = el[i]->GetDof();
+      order[i] = el[i]->GetOrder();
+      vdim[i] = elfun[i]->Size() / dof[i];
+
+      MFEM_ASSERT(vdim[i] == 1 ? true : hasFlag(modes_arr[i], ADEval::VECTOR),
+                  "ADNonlinearFormIntegrator: "
+                  "vdim must be 1 or the mode must be VECTOR");
+
+      elvect[i]->SetSize(elfun[i]->Size());
+      *elvect[i] = 0.0;
+   }
+
+   std::array<int, numSpaces> shapedim(InitInputShapes(el, Tr, allshapes));
+   Array<int> x_idx(numSpaces+1);
+   x_idx[0] = 0;
+   for (int i=0; i<numSpaces; i++)
+   {
+      x_idx[i+1] = shapedim[i]*vdim[i];
+   }
+   x_idx.PartialSum();
+   x.SetSize(f.n_input);
+   jac.SetSize(f.n_input);
+   _constexpr_for([&](auto vi)
+   {
+      xvar[vi].MakeRef(x, x_idx[vi], shapedim[vi]*vdim[vi]);
+      jacVar[vi].MakeRef(jac, x_idx[vi], shapedim[vi]*vdim[vi]);
+      if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+      {
+         elfun_matview[vi].UseExternalData(const_cast<real_t*>(elfun[vi]->GetData()),
+                                           dof[vi], vdim[vi]);
+         xmat[vi].UseExternalData(xvar[vi].GetData(), shapedim[vi], vdim[vi]);
+         jacVarMat[vi].UseExternalData(jacVar[vi].GetData(), shapedim[vi], vdim[vi]);
+      }
+   }, std::make_index_sequence<sizeof...(modes)> {});
+
+   const IntegrationRule * ir = GetIntegrationRule(el, Tr);
+   real_t w;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      w = Tr.Weight()*ip.weight;
+
+      CalcInputShapes(el, Tr, ip, allshapes);
+
+      _constexpr_for([&](auto vi)
+      {
+         if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+         {
+            MultAtB(allshapes[vi], elfun_matview[vi], xmat[vi]);
+         }
+         else
+         {
+            allshapes[vi].MultTranspose(*elfun[vi], xvar[vi]);
+         }
+      }, std::make_index_sequence<sizeof...(modes)> {});
+      f.Gradient(x, Tr, ip, jac);
+      jac *= w;
+
+      _constexpr_for([&](auto vi)
+      {
+         if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+         {
+            elvectmat[vi].UseExternalData(elvect[vi]->GetData(), dof[vi], vdim[vi]);
+            AddMult(allshapes[vi], jacVarMat[vi], elvectmat[vi]);
+         }
+         else
+         {
+            allshapes[vi].AddMult(jacVar[vi], *elvect[vi]);
+         }
+      }, std::make_index_sequence<sizeof...(modes)> {});
+   }
+}
+
+/// Perform the local action of the NonlinearFormIntegrator
+template <ADEval... modes>
+void ADBlockNonlinearFormIntegrator<modes...>::AssembleElementGrad(
+   const Array<const FiniteElement *>&el,
+   ElementTransformation &Tr,
+   const Array<const Vector *>&elfun,
+   const Array2D<DenseMatrix *>&elmat)
+{
+   MFEM_ASSERT(el.Size() == numSpaces,
+               "ADBlockNonlinearFormIntegrator: "
+               "el.Size()=" << el.Size() << " must match numSpaces=" << numSpaces);
+   Array<int> dof(numSpaces);
+   Array<int> order(numSpaces);
+   for (int i=0; i<numSpaces; i++)
+   {
+      dof[i] = el[i]->GetDof();
+      order[i] = el[i]->GetOrder();
+      vdim[i] = elfun[i]->Size() / dof[i];
+
+      MFEM_ASSERT(vdim[i] == 1 ? true : hasFlag(modes_arr[i], ADEval::VECTOR),
+                  "ADNonlinearFormIntegrator: "
+                  "vdim must be 1 or the mode must be VECTOR");
+   }
+
+   for (int j=0; j<numSpaces; j++)
+   {
+      for (int i=0; i<numSpaces; i++)
+      {
+         elmat(i,j)->SetSize(elfun[i]->Size(),
+                             elfun[j]->Size());
+         *elmat(i,j) = 0.0;
+      }
+   }
+
+   std::array<int, numSpaces> shapedim(InitInputShapes(el, Tr, allshapes));
+   Array<int> x_idx(numSpaces+1);
+   x_idx[0] = 0;
+   for (int i=0; i<numSpaces; i++)
+   {
+      x_idx[i+1] = shapedim[i]*vdim[i];
+   }
+   x_idx.PartialSum();
+   x.SetSize(f.n_input);
+   H.SetSize(f.n_input);
+
+   _constexpr_for([&](auto vi)
+   {
+      xvar[vi].MakeRef(x, x_idx[vi], shapedim[vi]*vdim[vi]);
+      if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+      {
+         elfun_matview[vi].UseExternalData(const_cast<real_t*>(elfun[vi]->GetData()),
+                                           dof[vi], vdim[vi]);
+         xmat[vi].UseExternalData(xvar[vi].GetData(), shapedim[vi], vdim[vi]);
+      }
+   }, std::make_index_sequence<sizeof...(modes)> {});
+
+   const IntegrationRule * ir = GetIntegrationRule(el, Tr);
+   real_t w;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      w = Tr.Weight()*ip.weight;
+
+      CalcInputShapes(el, Tr, ip, allshapes);
+
+      _constexpr_for([&](auto vi)
+      {
+         if constexpr (hasFlag(modes_arr[vi], ADEval::VECTOR))
+         {
+            MultAtB(allshapes[vi], elfun_matview[vi], xmat[vi]);
+         }
+         else
+         {
+            allshapes[vi].MultTranspose(*elfun[vi], xvar[vi]);
+         }
+      }, std::make_index_sequence<sizeof...(modes)> {});
+      f.Hessian(x, Tr, ip, H);
+      H *= w;
+      _constexpr_for([&](auto trial_i)
+      {
+         _constexpr_for([&](auto test_i)
+         {
+            const int tr_vdim = vdim[trial_i];
+            const int ts_vdim = vdim[test_i];
+            H.GetSubMatrix(x_idx[test_i], x_idx[test_i+1], x_idx[trial_i], x_idx[trial_i+1],
+                           Hsub);
+            Hsub.SetSize(shapedim[test_i], ts_vdim*tr_vdim*shapedim[trial_i]);
+            Hx.SetSize(dof[test_i], ts_vdim*tr_vdim*shapedim[trial_i]);
+            Mult(allshapes[test_i], Hsub, Hx);
+            Hx.SetSize(dof[test_i]*ts_vdim, tr_vdim*shapedim[trial_i]);
+            const int h = dof[test_i]*ts_vdim;
+            const int w = shapedim[trial_i];
+            const int wout = dof[trial_i];
+            for (int d=0; d<tr_vdim; d++)
+            {
+               Hxsub.UseExternalData(Hx.GetData() + d*(w*h), h, w);
+               partelmat[trial_i].UseExternalData(elmat(test_i, trial_i)->GetData() + d*wout*h,
+                                                  h, wout);
+               AddMultABt(Hxsub, allshapes[trial_i], partelmat[trial_i]);
+            }
+         }, std::make_index_sequence<sizeof...(modes)> {});
+      }, std::make_index_sequence<sizeof...(modes)> {});
+   }
+}
+
+/// @brief Perform the local action of the NonlinearFormIntegrator resulting
+/// from a face integral term.
+template <ADEval... modes>
+void ADBlockNonlinearFormIntegrator<modes...>::AssembleFaceVector(
+   const Array<const FiniteElement *>&el1,
+   const Array<const FiniteElement *>&el2,
+   FaceElementTransformations &Tr,
+   const Array<const Vector *>&elfun,
+   const Array<Vector *>&elvect)
+{
+   MFEM_ABORT("ADBlockNonlinearFormIntegrator::AssembleFaceVector: "
+              "This method is not implemented.");
+}
+
+
+/// @brief Assemble the local action of the gradient of the
+/// NonlinearFormIntegrator resulting from a face integral term.
+template <ADEval... modes>
+void ADBlockNonlinearFormIntegrator<modes...>::AssembleFaceGrad(
+   const Array<const FiniteElement *>&el1,
+   const Array<const FiniteElement *>&el2,
+   FaceElementTransformations &Tr,
+   const Array<const Vector *>&elfun,
+   const Array2D<DenseMatrix *>&elmat)
+{
+   MFEM_ABORT("ADBlockNonlinearFormIntegrator::AssembleFaceGrad: "
+              "This method is not implemented.");
+}
+} // namespace mfem

miniapps/autodiff/ad_ex0.cpp

@@ -0,0 +1,187 @@
+/// Example 0: AD Function Example
+#include "mfem.hpp"
+#include "ad_native.hpp"
+
+using namespace std;
+using namespace mfem;
+
+void Print(DenseMatrix &H, bool last_linebreak=true);
+void Print(DenseTensor &H);
+
+struct MyADFunction : public ADFunction
+{
+public:
+   MyADFunction(int n_input): ADFunction(n_input) { }
+
+   AD_IMPL(T, V, M, x, return sin(x(0))*exp(x(1)) + pow(x(2), 3.0);)
+};
+
+struct MyADVecFunction : public ADVectorFunction
+{
+public:
+   MyADVecFunction()
+      : ADVectorFunction(3, 2) { }
+
+   AD_VEC_IMPL(T, V, M, x, result,
+   {
+      result[0]=sin(x[0]*x[1]);
+      result[1]=cos(x[0]*x[1]*x[2]);
+   });
+};
+
+void jacobian(const Vector &x, Vector &J)
+{
+   J.SetSize(x.Size());
+   J[0] = std::cos(x(0)) * std::exp(x(1));
+   J[1] = std::sin(x(0)) * std::exp(x(1));
+   J[2] = 3.0 * std::pow(x(2), 2.0);
+}
+
+void hessian(const Vector &x, DenseMatrix &H)
+{
+   // J[0] = cos(x(0)) * exp(x(1));
+   H.SetSize(x.Size(), x.Size());
+   H(0, 0) = -std::sin(x(0)) * std::exp(x(1));
+   H(0, 1) = std::cos(x(0)) *std::exp(x(1));
+   H(0, 2) = 0.0;
+
+   // J[1] = sin(x(0)) * exp(x(1));
+   H(1, 0) = std::cos(x(0)) * std::exp(x(1));
+   H(1, 1) = std::sin(x(0)) * std::exp(x(1));
+   H(1, 2) = 0.0;
+
+   // J[2] = 3.0 * pow(x(2), 2.0);
+   H(2, 0) = 0.0;
+   H(2, 1) = 0.0;
+   H(2, 2) = 6.0 * std::pow(x(2), 1.0);
+}
+
+void jacobian(const Vector &x, DenseMatrix &J)
+{
+   // result[0]=sin(x[0]*x[1]);
+   // result[1]=cos(x[0]*x[1]*x[2]);)
+   J.SetSize(2,3);
+   J(0,0) = x(1) * std::cos(x(0) * x(1));
+   J(0,1) = x(0) * std::cos(x(0) * x(1));
+   J(0,2) = 0.0;
+   J(1,0) = -x(1) * x(2) * std::sin(x(0) * x(1) * x(2));
+   J(1,1) = -x(0) * x(2) * std::sin(x(0) * x(1) * x(2));
+   J(1,2) = -x(0) * x(1) * std::sin(x(0) * x(1) * x(2));
+}
+void hessian(const Vector &X, DenseTensor &H)
+{
+   real_t x(X(0)), y(X(1)), z(X(2));
+   H.SetSize(3, 3, 2);
+   H = 0.0;
+   using std::sin;
+   using std::cos;
+
+   // result[0]=sin(x[0]*x[1]);
+   H(0,0,0) = -y*y*sin(x*y);
+   H(0,1,0) = cos(x*y) - x*y*sin(x*y);
+   H(1,0,0) = cos(x*y) - x*y*sin(x*y);
+   H(1,1,0) = -x*x*sin(x*y);
+   // result[1]=cos(x[0]*x[1]*x[2]);)
+   H(0,0,1) = -y*y*z*z*cos(x*y*z);
+   H(1,0,1) = -x*y*z*z*cos(x*y*z) - z*sin(x*y*z);
+   H(2,0,1) = -x*y*y*z*cos(x*y*z) - y*sin(x*y*z);
+   H(0,1,1) = -x*y*z*z*cos(x*y*z) - z*sin(x*y*z);
+   H(1,1,1) = -x*x*z*z*cos(x*y*z);
+   H(2,1,1) = -x*x*y*z*cos(x*y*z) - x*sin(x*y*z);
+   H(0,2,1) = -x*y*y*z*cos(x*y*z) - y*sin(x*y*z);
+   H(1,2,1) = -x*x*y*z*cos(x*y*z) - x*sin(x*y*z);
+   H(2,2,1) = -x*x*y*y*cos(x*y*z);
+}
+
+int main(int argc, char *argv[])
+{
+   Vector x({0.5, 1.0, -1.0});
+
+   MyADFunction f(3);
+
+   Vector jac, jac_ref;
+   f.Gradient(x, jac);
+   jacobian(x, jac_ref);
+
+   DenseMatrix hess, hess_ref;
+   f.Hessian(x, hess);
+   hessian(x, hess_ref);
+
+   MyADVecFunction f2;
+
+   DenseMatrix jac2, jac2_ref;
+   f2.Gradient(x, jac2);
+   jacobian(x, jac2_ref);
+
+   DenseTensor hess2, hess2_ref;
+   f2.Hessian(x, hess2);
+   hessian(x, hess2_ref);
+
+
+   out << "Value : " << f(x) << std::endl;
+
+   out << "Jacobian  : ";
+   jac.Print();
+   out << "Reference : ";
+   jac_ref.Print();
+   jac -= jac_ref;
+
+   out << "Hessian : " << std::endl;
+   Print(hess);
+   out << "Reference: " << std::endl;
+   Print(hess_ref);
+   hess -= hess_ref;
+
+   out << std::endl;
+   out << "Jacobian error: " << jac.DistanceTo(jac_ref) << std::endl;
+   out << "Hessian error: " << hess.MaxMaxNorm() << std::endl;
+   out << "-------------------------" << std::endl;
+
+   out << "Jacobian2 : " << std::endl;
+   Print(jac2);
+   out << "Reference : " << std::endl;
+   Print(jac2_ref);
+   jac2 -= jac2_ref;
+
+   out << "Hess2 : " << std::endl;
+   Print(hess2);
+   out << "Reference : " << std::endl;
+   Print(hess2_ref);
+   out << std::endl;
+   out << "Jacobian2 error: " << jac2.MaxMaxNorm() << std::endl;
+   for (int k=0; k<hess2.SizeK(); k++)
+   {
+      hess2(k) -= hess2_ref(k);
+      out << "Hessian[" << k << "] error: " << hess2(k).MaxMaxNorm() << std::endl;
+   }
+
+   return 0;
+}
+
+void Print(DenseMatrix &H, bool last_linebreak)
+{
+   for (int i = 0; i < H.Height(); i++)
+   {
+      for (int j = 0; j < H.Width(); j++)
+      {
+         cout << H(i, j) << " ";
+      }
+      cout << ";";
+      if (i < H.Height() - 1 || last_linebreak)
+      {
+         cout << "\n";
+      }
+   }
+   cout << std::flush;
+}
+
+void Print(DenseTensor &H)
+{
+   for (int k=0; k<H.SizeK(); k++)
+   {
+      out << "{ ";
+      Print(H(k), false);
+      out << " }\n";
+   }
+   cout << std::flush;
+}

miniapps/autodiff/ad_ex1.cpp

@@ -0,0 +1,96 @@
+/// Example 1: AD Diffusion
+#include "mfem.hpp"
+#include "logger.hpp"
+#include "ad_intg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // file name to be saved
+   std::stringstream filename;
+   filename << "ad-diffusion";
+
+   int order = 1;
+   int ref_levels = 1;
+   bool visualization = false;
+   bool paraview = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ref_levels, "-r", "--ref", "Refinement levels");
+   args.AddOption(&visualization, "-vis", "--visualization",
+                  "-no-vis", "--no-visualization",
+                  "Enable visualization, default is false");
+   args.AddOption(&paraview, "-pv", "--paraview",
+                  "-no-pv", "--no-paraview",
+                  "Enable Paraview Export. Default is false");
+   args.ParseCheck();
+
+   Mesh mesh = Mesh::MakeCartesian2D(10, 10,
+                                     Element::QUADRILATERAL);
+   const int dim = mesh.Dimension();
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+   FunctionCoefficient load_cf([](const Vector &x)
+   {
+      return 2*M_PI * M_PI * std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
+   });
+
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   Array<int> ess_tdof_list;
+   fes.GetBoundaryTrueDofs(ess_tdof_list);
+
+   DiffusionEnergy energy(dim);
+
+   NonlinearForm nlf(&fes);
+   nlf.AddDomainIntegrator(new ADNonlinearFormIntegrator<ADEval::GRAD>(energy));
+   nlf.SetEssentialTrueDofs(ess_tdof_list);
+   LinearForm load(&fes);
+   load.AddDomainIntegrator(new DomainLFIntegrator(load_cf));
+   load.Assemble();
+   load.SetSubVector(ess_tdof_list, 0.0);
+
+   GridFunction x(&fes);
+   x = 0.0;
+   SparseMatrix &op = static_cast<SparseMatrix&>(nlf.GetGradient(x));
+   CGSolver lin_solver;
+   GSSmoother prec;
+   lin_solver.SetPreconditioner(prec);
+   lin_solver.SetOperator(op);
+   lin_solver.SetRelTol(1e-12);
+   lin_solver.SetAbsTol(0.0);
+   lin_solver.SetMaxIter(1e04);
+   lin_solver.Mult(load, x);
+   if (visualization)
+   {
+      GLVis glvis("localhost", 19916);
+      glvis.Append(x, "x", "Rjc");
+   }
+   if (paraview)
+   {
+      std::stringstream pvloc;
+      pvloc << "ParaView/" << filename.str();
+      ParaViewDataCollection paraview_dc(pvloc.str(), &mesh);
+      paraview_dc.SetLevelsOfDetail(order);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.RegisterField("solution", &x);
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetTime(0.0);
+      paraview_dc.Save();
+   }
+
+   FunctionCoefficient exact_sol([](const Vector &x)
+   {
+      return std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
+   });
+   real_t err = x.ComputeL2Error(exact_sol);
+   out << "Error: " << err << std::endl;
+   return 0;
+}

miniapps/autodiff/ad_ex2.cpp

@@ -0,0 +1,119 @@
+/// Example 2: AD Minimal Surface
+#include "mfem.hpp"
+#include "logger.hpp"
+#include "ad_intg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+struct MinimalSurfaceEnergy : public ADFunction
+{
+public:
+   real_t eps=0.5; // regularization
+   MinimalSurfaceEnergy(int dim): ADFunction(dim) {}
+   AD_IMPL(T, V, M, gradu,
+   {
+      T h1_norm(gradu*gradu);
+      // sqrt(1+ ||grad u||^2)
+      // dJ/du = 0 -> minimal surface
+      return sqrt(h1_norm + 1.0) + eps*h1_norm;
+   });
+};
+
+int main(int argc, char *argv[])
+{
+   // file name to be saved
+   std::stringstream filename;
+   filename << "ad-minimalsurface";
+
+   int order = 1;
+   int ref_levels = 3;
+   bool visualization = false;
+   bool paraview = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ref_levels, "-r", "--ref", "Refinement levels");
+   args.AddOption(&visualization, "-vis", "--visualization",
+                  "-no-vis", "--no-visualization",
+                  "Enable visualization, default is false");
+   args.AddOption(&paraview, "-pv", "--paraview",
+                  "-no-pv", "--no-paraview",
+                  "Enable Paraview Export. Default is false");
+   args.ParseCheck();
+
+   // Mesh mesh = rhs_fun_circle
+   Mesh mesh = Mesh::MakeCartesian2D(10, 10,
+                                     Element::QUADRILATERAL);
+   const int dim = mesh.Dimension();
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+   FunctionCoefficient bdry_cf([](const Vector &x)
+   {
+      real_t theta = std::atan2(x(1)-0.5, x(0)-0.5);
+      real_t r = std::sqrt(std::pow(x(0)-0.5, 2.0) + std::pow(x(1)-0.5, 2.0));
+      return r*std::cos(2*theta);
+   });
+
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   Array<int> is_bdr_ess(mesh.bdr_attributes.Max());
+   is_bdr_ess = 1;
+
+   MinimalSurfaceEnergy energy(dim);
+
+   NonlinearForm nlf(&fes);
+   nlf.AddDomainIntegrator(new ADNonlinearFormIntegrator<ADEval::GRAD>
+                           (energy));
+   nlf.SetEssentialBC(is_bdr_ess);
+
+   GridFunction x(&fes);
+   x = 0.0;
+   x.ProjectBdrCoefficient(bdry_cf, is_bdr_ess);
+   CGSolver lin_solver;
+   GSSmoother prec;
+   lin_solver.SetPreconditioner(prec);
+   lin_solver.SetRelTol(1e-12);
+   lin_solver.SetAbsTol(0.0);
+   lin_solver.SetMaxIter(1e04);
+   NewtonSolver solver;
+   solver.SetSolver(lin_solver);
+   solver.SetOperator(nlf);
+   solver.SetAbsTol(1e-10);
+   solver.SetRelTol(1e-10);
+   IterativeSolver::PrintLevel print_level;
+   print_level.iterations = 1;
+   solver.SetPrintLevel(print_level);
+   solver.SetMaxIter(100);
+   solver.iterative_mode = true;
+   Vector dummy(0);
+   std::unique_ptr<GLVis> glvis;
+   if (visualization)
+   {
+      glvis = std::make_unique<GLVis>("localhost", 19916);
+      glvis->Append(x, "x", "Rjc");
+   }
+   for (int i=0; i<30; i++)
+   {
+      solver.Mult(dummy, x);
+      if (glvis) { glvis->Update(); }
+      energy.eps *= 0.5;
+   }
+   if (paraview)
+   {
+      std::stringstream pvloc;
+      pvloc << "ParaView/" << filename.str();
+      ParaViewDataCollection paraview_dc(pvloc.str(), &mesh);
+      paraview_dc.SetLevelsOfDetail(order);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.RegisterField("solution", &x);
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetTime(0.0);
+      paraview_dc.Save();
+   }
+   return 0;
+}

miniapps/autodiff/ad_ex3.cpp

@@ -0,0 +1,98 @@
+/// Example 3: AD Linear Elasticity with Vector FE
+#include "mfem.hpp"
+#include "logger.hpp"
+#include "ad_intg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   // file name to be saved
+   std::stringstream filename;
+   filename << "ad-elasticity";
+
+   int order = 1;
+   int ref_levels = 3;
+   bool visualization = false;
+   bool paraview = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ref_levels, "-r", "--ref", "Refinement levels");
+   args.AddOption(&visualization, "-vis", "--visualization",
+                  "-no-vis", "--no-visualization",
+                  "Enable visualization, default is false");
+   args.AddOption(&paraview, "-pv", "--paraview",
+                  "-no-pv", "--no-paraview",
+                  "Enable Paraview Export. Default is false");
+   args.ParseCheck();
+
+   // Mesh mesh = rhs_fun_circle
+   Mesh mesh = Mesh::MakeCartesian2D(10, 10,
+                                     Element::QUADRILATERAL);
+   const int dim = mesh.Dimension();
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+   VectorFunctionCoefficient load_cf(dim, [dim](const Vector &x, Vector &y)
+   {
+      y.SetSize(dim);
+      y = 1.0;
+   });
+
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec, dim);
+   Array<int> is_bdr_ess(mesh.bdr_attributes.Max());
+   is_bdr_ess = 0;
+   is_bdr_ess[3] = 1;
+   Array<int> ess_tdof_list;
+   fes.GetEssentialTrueDofs(is_bdr_ess, ess_tdof_list);
+
+   real_t lambda(1.0), mu(1.0);
+   LinearElasticityEnergy energy(dim, lambda, mu);
+
+   NonlinearForm nlf(&fes);
+   nlf.AddDomainIntegrator(
+      new ADNonlinearFormIntegrator<ADEval::GRAD | ADEval::VECTOR>(energy));
+   nlf.SetEssentialBC(is_bdr_ess);
+   LinearForm load(&fes);
+   load.AddDomainIntegrator(new VectorDomainLFIntegrator(load_cf));
+   load.Assemble();
+   load.SetSubVector(ess_tdof_list, 0.0);
+
+   GridFunction x(&fes);
+   x = 0.0;
+   SparseMatrix &op = static_cast<SparseMatrix&>(nlf.GetGradient(x));
+   CGSolver lin_solver;
+   GSSmoother prec;
+   lin_solver.SetPreconditioner(prec);
+   lin_solver.SetOperator(op);
+   lin_solver.SetRelTol(1e-12);
+   lin_solver.SetAbsTol(0.0);
+   lin_solver.SetMaxIter(1e04);
+   lin_solver.Mult(load, x);
+
+   if (visualization)
+   {
+      GLVis glvis("localhost", 19916);
+      glvis.Append(x, "x", "Rjc");
+   }
+   if (paraview)
+   {
+      std::stringstream pvloc;
+      pvloc << "ParaView/" << filename.str();
+      ParaViewDataCollection paraview_dc(pvloc.str(), &mesh);
+      paraview_dc.SetLevelsOfDetail(order);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.RegisterField("solution", &x);
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetTime(0.0);
+      paraview_dc.Save();
+   }
+   return 0;
+}

miniapps/autodiff/ad_ex4.cpp

@@ -0,0 +1,232 @@
+/// Example 4: AD Obstacle Problem with PG
+#include "mfem.hpp"
+#include "logger.hpp"
+#include "ad_intg.hpp"
+#include "tools.hpp"
+#include "pg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+struct ObstacleEnergy : public ADFunction
+{
+   ObstacleEnergy(int dim) : ADFunction(dim+1) {}
+   AD_IMPL(T, V, M, x,
+   {
+      T result = {};
+      // First component is u. Others are grad u
+      for (int i=1; i<x.Size(); i++)
+      {
+         result += x[i]*x[i];
+      }
+      return result*0.5;
+   });
+};
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+   MPI_Comm comm = MPI_COMM_WORLD;
+   // file name to be saved
+   std::stringstream filename;
+   filename << "ad-obstacle";
+   int rule_type = PGStepSizeRule::RuleType::CONSTANT;
+   real_t max_alpha = 1e04;
+   real_t alpha0 = 1.0;
+   real_t ratio = 1.0;
+   real_t ratio2 = 1.0;
+
+   int order = 2;
+   int ref_levels = 3;
+   bool visualization = false;
+   bool paraview = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ref_levels, "-r", "--ref", "Refinement levels");
+   args.AddOption(&rule_type, "-rule", "--rule",
+                  "Step size rule type: 0=CONSTANT, 1=POLY, 2=EXP, 3=DOUBLE_EXP");
+   args.AddOption(&max_alpha, "-ma", "--max-alpha",
+                  "Maximum step size for PG method");
+   args.AddOption(&alpha0, "-a0", "--alpha0",
+                  "Initial step size for PG method");
+   args.AddOption(&ratio, "-ar", "--alpha-ratio",
+                  "Ratio for step size rule (POLY, EXP, DOUBLE_EXP)");
+   args.AddOption(&ratio2, "-ar2", "--alpha-ratio2",
+                  "Second ratio for DOUBLE_EXP step size rule");
+   args.AddOption(&visualization, "-vis", "--visualization",
+                  "-no-vis", "--no-visualization",
+                  "Enable visualization, default is false");
+   args.AddOption(&paraview, "-pv", "--paraview",
+                  "-no-pv", "--no-paraview",
+                  "Enable Paraview Export. Default is false");
+   args.ParseCheck();
+   if (myid != 0) { out.Disable(); }
+
+   PGStepSizeRule alpha_rule(rule_type, alpha0, max_alpha, ratio, ratio2);
+
+   // Mesh mesh = rhs_fun_circle
+   Mesh ser_mesh = Mesh::MakeCartesian2D(2, 2,
+                                         Element::QUADRILATERAL);
+   const int dim = ser_mesh.Dimension();
+   for (int i = 0; i < ref_levels; i++)
+   {
+      ser_mesh.UniformRefinement();
+   }
+   ParMesh mesh(comm, ser_mesh);
+
+   const int numBdrAttr = mesh.bdr_attributes.Max();
+   Array<int> is_bdr_ess1(numBdrAttr);
+   is_bdr_ess1 = 1;
+   Array<int> is_bdr_ess2(numBdrAttr);
+   is_bdr_ess2 = 0;
+   Array<Array<int>*> is_bdr_ess{&is_bdr_ess1, &is_bdr_ess2};
+   FunctionCoefficient load_cf([](const Vector &x)
+   {
+      return 2*M_PI * M_PI * std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
+   });
+   ObstacleEnergy obj_energy(dim);
+
+   H1_FECollection primal_fec(order+1, dim);
+   L2_FECollection latent_fec(order-1, dim);
+   ParFiniteElementSpace primal_fes(&mesh, &primal_fec);
+   ParFiniteElementSpace latent_fes(&mesh, &latent_fec);
+   QuadratureSpace visspace(&mesh, order+3);
+   const IntegrationRule &ir = IntRules.Get(Geometry::Type::SQUARE, 3*order + 3);
+
+   Array<int> ess_tdof_list;
+   primal_fes.GetEssentialTrueDofs(is_bdr_ess1, ess_tdof_list);
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = primal_fes.GetTrueVSize();
+   offsets[2] = latent_fes.GetTrueVSize();
+   offsets.PartialSum();
+   BlockVector x_and_latent(offsets);
+
+   ParGridFunction x(&primal_fes), latent(&latent_fes);
+   ParGridFunction latent_k(latent);
+
+   x = 0.0; x.ParallelAssemble(x_and_latent.GetBlock(0));
+   latent = 0.0; latent.ParallelAssemble(x_and_latent.GetBlock(1));
+   latent_k = 0.0; latent_k.SetTrueVector();
+
+   FermiDiracEntropy entropy(0.0, 0.5);
+
+   DifferentiableCoefficient entropy_cf(entropy);
+   entropy_cf.AddInput(&latent);
+   VectorCoefficient &u_cf = entropy_cf.Gradient();
+
+   real_t alpha;
+   ADPGFunctional pg_functional(obj_energy, entropy, &alpha, latent_k);
+
+   ParGridFunction lambda(latent), lambda_prev(latent);
+   lambda = 0.0;
+   GridFunctionCoefficient lambda_prev_cf(&lambda_prev);
+
+   Array<ParFiniteElementSpace*> fespaces{&primal_fes, &latent_fes};
+   ParBlockNonlinearForm bnlf(fespaces);
+   constexpr ADEval u_mode = ADEval::VALUE | ADEval::GRAD;
+   constexpr ADEval latent_mode = ADEval::VALUE;
+   bnlf.AddDomainIntegrator(
+      new ADBlockNonlinearFormIntegrator<u_mode, latent_mode>(
+         pg_functional, &ir)
+   );
+
+   BlockVector rhs(offsets);
+   ParLinearForm b(&primal_fes);
+   b.AddDomainIntegrator(new DomainLFIntegrator(load_cf));
+   b.Assemble();
+   b.ParallelAssemble(rhs.GetBlock(0));
+   rhs.GetBlock(0).SetSubVector(ess_tdof_list, 0.0);
+   rhs.GetBlock(1) = 0.0;
+
+   Array<Vector*> rhs_list{&rhs.GetBlock(0), &rhs.GetBlock(1)};
+   bnlf.SetEssentialBC(is_bdr_ess, rhs_list);
+
+   MUMPSMonoSolver lin_solver(comm);
+   NewtonSolver solver(comm);
+   solver.SetSolver(lin_solver);
+   solver.SetOperator(bnlf);
+   IterativeSolver::PrintLevel print_level;
+   solver.SetPrintLevel(print_level);
+   solver.SetAbsTol(1e-09);
+   solver.SetRelTol(0.0);
+   solver.SetMaxIter(20);
+   solver.iterative_mode = true;
+
+   std::unique_ptr<GLVis> glvis;
+   if (visualization)
+   {
+      glvis = std::make_unique<GLVis>("localhost", 19916, 400, 350, 3);
+      glvis->Append(x, "u", "Rjclmm");
+      glvis->Append(u_cf, visspace, "U(psi)", "RjclQmm");
+      glvis->Append(lambda, "lambda", "Rjclmm");
+   }
+   std::unique_ptr<ParaViewDataCollection> paraview_dc;
+   if (paraview)
+   {
+      filename << "r" << ref_levels << "-o" << order;
+      paraview_dc = std::make_unique<ParaViewDataCollection>(filename.str(), &mesh);
+      paraview_dc->SetLevelsOfDetail(order);
+      paraview_dc->SetDataFormat(VTKFormat::BINARY);
+      paraview_dc->SetHighOrderOutput(true);
+      paraview_dc->RegisterField("solution", &x);
+      paraview_dc->SetCycle(0);
+      paraview_dc->SetTime(0.0);
+      paraview_dc->Save();
+   }
+
+
+   real_t lambda_diff = infinity();
+   for (int i=0; i<100; i++)
+   {
+      alpha = alpha_rule.Get(i);
+      out << "PG iteration " << i + 1 << " with alpha=" << alpha << std::endl;
+      latent_k = latent;
+      latent_k.SetTrueVector();
+
+      solver.Mult(rhs, x_and_latent);
+
+      if (!solver.GetConverged())
+      {
+         out << "Newton Failed to converge in " << solver.GetNumIterations() <<
+             std::endl;
+      }
+      x.SetFromTrueDofs(x_and_latent.GetBlock(0));
+      latent.SetFromTrueDofs(x_and_latent.GetBlock(1));
+
+      if (glvis) { glvis->Update(); }
+      if (paraview_dc)
+      {
+         paraview_dc->SetCycle(i+1);
+         paraview_dc->SetTime(i+1);
+         paraview_dc->Save();
+      }
+
+      subtract(latent, latent_k, lambda);
+      lambda *= 1.0 / pg_functional.GetAlpha();
+
+      if ((lambda_diff = lambda.ComputeL1Error(lambda_prev_cf)) < 1e-8)
+      {
+         out << "  The dual variable, (psi - psi_k)/alpha, converged" << std::endl;
+         out << "PG Converged in " << i + 1
+             << " with final Lambda difference: " << lambda_diff << std::endl;
+         break;
+      }
+      else
+      {
+         out << "  Newton converged in " << solver.GetNumIterations()
+             << " with residual " << solver.GetFinalNorm() << std::endl;
+         out << "  Lambda difference: " << lambda_diff << std::endl;
+      }
+
+      lambda_prev = lambda;
+   }
+   return 0;
+}

miniapps/autodiff/ad_ex5.cpp

@@ -0,0 +1,232 @@
+/// Example 5: AD Gradeint Obstacle Problem with PG
+#include "mfem.hpp"
+#include "logger.hpp"
+#include "ad_intg.hpp"
+#include "tools.hpp"
+#include "pg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+struct ObstacleEnergy : public ADFunction
+{
+   ObstacleEnergy(int dim) : ADFunction(dim) {}
+   AD_IMPL(T, V, M, x,
+   {
+      return x*x*0.5;
+   });
+};
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+   MPI_Comm comm = MPI_COMM_WORLD;
+   // file name to be saved
+   std::stringstream filename;
+   filename << "ad-grad-obstacle";
+   int rule_type = PGStepSizeRule::RuleType::CONSTANT;
+   real_t max_alpha = 1e06;
+   real_t alpha0 = 1.0;
+   real_t ratio = 1.0;
+   real_t ratio2 = 1.0;
+   bool use_iterative = false;
+
+   int order = 2;
+   int ref_levels = 3;
+   bool visualization = false;
+   bool paraview = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ref_levels, "-r", "--ref", "Refinement levels");
+   args.AddOption(&rule_type, "-rule", "--rule",
+                  "Step size rule type: 0=CONSTANT, 1=POLY, 2=EXP, 3=DOUBLE_EXP");
+   args.AddOption(&max_alpha, "-ma", "--max-alpha",
+                  "Maximum step size for PG method");
+   args.AddOption(&alpha0, "-a0", "--alpha0",
+                  "Initial step size for PG method");
+   args.AddOption(&ratio, "-ar", "--alpha-ratio",
+                  "Ratio for step size rule (POLY, EXP, DOUBLE_EXP)");
+   args.AddOption(&ratio2, "-ar2", "--alpha-ratio2",
+                  "Second ratio for DOUBLE_EXP step size rule");
+   args.AddOption(&visualization, "-vis", "--visualization",
+                  "-no-vis", "--no-visualization",
+                  "Enable visualization, default is false");
+   args.AddOption(&paraview, "-pv", "--paraview",
+                  "-no-pv", "--no-paraview",
+                  "Enable Paraview Export. Default is false");
+   args.AddOption(&use_iterative, "-gmres", "--preconditioned-gmres",
+                  "-mumps", "--MUMPS",
+                  "Use preconditioned GMRES or MUMPS as linear solver. Default is MUMPS");
+   args.ParseCheck();
+   if (myid != 0) { out.Disable(); }
+   MFEMInitializePetsc(NULL,NULL,"../src/pgpetsc",NULL);
+
+   PGStepSizeRule alpha_rule(rule_type, alpha0, max_alpha, ratio, ratio2);
+
+   // Mesh mesh = rhs_fun_circle
+   Mesh ser_mesh = Mesh::MakeCartesian2D(2, 2,
+                                         Element::QUADRILATERAL);
+   const int dim = ser_mesh.Dimension();
+   for (int i = 0; i < ref_levels; i++)
+   {
+      ser_mesh.UniformRefinement();
+   }
+   ParMesh mesh(comm, ser_mesh);
+
+   const int numBdrAttr = mesh.bdr_attributes.Max();
+   Array<int> is_bdr_ess1(numBdrAttr);
+   is_bdr_ess1 = 1;
+   Array<int> is_bdr_ess2(numBdrAttr);
+   is_bdr_ess2 = 0;
+   Array<Array<int>*> is_bdr_ess{&is_bdr_ess1, &is_bdr_ess2};
+   FunctionCoefficient load_cf([](const Vector &x)
+   {
+      return 2*M_PI * M_PI * std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
+   });
+   ObstacleEnergy obj_energy(dim);
+
+   H1_FECollection primal_fec(order, dim);
+   H1_FECollection latent_fec(order-1, dim);
+   ParFiniteElementSpace primal_fes(&mesh, &primal_fec);
+   ParFiniteElementSpace latent_fes(&mesh, &latent_fec, dim);
+   QuadratureSpace visspace(&mesh, order+3);
+   const IntegrationRule &ir = IntRules.Get(Geometry::Type::SQUARE, 3*order + 3);
+
+   Array<int> ess_tdof_list;
+   primal_fes.GetEssentialTrueDofs(is_bdr_ess1, ess_tdof_list);
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = primal_fes.GetTrueVSize();
+   offsets[2] = latent_fes.GetTrueVSize();
+   offsets.PartialSum();
+   BlockVector x_and_latent(offsets);
+
+   ParGridFunction x(&primal_fes), latent(&latent_fes);
+   ParGridFunction latent_k(latent);
+
+   x = 0.0; x.ParallelAssemble(x_and_latent.GetBlock(0));
+   latent = 0.0; latent.ParallelAssemble(x_and_latent.GetBlock(1));
+   latent_k = 0.0; latent_k.SetTrueVector();
+
+   FunctionCoefficient bound([](const Vector &x)
+   { return 0.1 + 0.2*x[0] + 0.4*x[1]; });
+   HellingerEntropy entropy(dim, &bound);
+
+   DifferentiableCoefficient entropy_cf(entropy);
+   entropy_cf.AddInput(&latent);
+   VectorCoefficient &u_cf = entropy_cf.Gradient();
+
+   real_t alpha;
+   ADPGFunctional pg_functional(obj_energy, entropy, &alpha, latent_k);
+
+   ParGridFunction lambda(latent), lambda_prev(latent);
+   lambda = 0.0;
+   VectorGridFunctionCoefficient lambda_prev_cf(&lambda_prev);
+
+   Array<ParFiniteElementSpace*> fespaces{&primal_fes, &latent_fes};
+   ParBlockNonlinearForm bnlf(fespaces);
+   constexpr ADEval u_mode = ADEval::GRAD;
+   constexpr ADEval latent_mode = ADEval::VALUE | ADEval::VECTOR;
+   bnlf.AddDomainIntegrator(
+      new ADBlockNonlinearFormIntegrator<u_mode, latent_mode>(
+         pg_functional, &ir)
+   );
+
+   BlockVector rhs(offsets);
+   ParLinearForm b(&primal_fes);
+   b.AddDomainIntegrator(new DomainLFIntegrator(load_cf));
+   b.Assemble();
+   b.ParallelAssemble(rhs.GetBlock(0));
+   rhs.GetBlock(0).SetSubVector(ess_tdof_list, 0.0);
+   rhs.GetBlock(1) = 0.0;
+
+   Array<Vector*> rhs_list{&rhs.GetBlock(0), &rhs.GetBlock(1)};
+   bnlf.SetEssentialBC(is_bdr_ess, rhs_list);
+
+   MUMPSMonoSolver lin_solver(comm);
+   NewtonSolver solver(comm);
+   solver.SetSolver(lin_solver);
+   solver.SetOperator(bnlf);
+   IterativeSolver::PrintLevel print_level;
+   solver.SetPrintLevel(print_level);
+   solver.SetAbsTol(1e-09);
+   solver.SetRelTol(0.0);
+   solver.SetMaxIter(20);
+   solver.iterative_mode = true;
+
+   std::unique_ptr<GLVis> glvis;
+   if (visualization)
+   {
+      glvis = std::make_unique<GLVis>("localhost", 19916, 400, 350, 3);
+      glvis->Append(x, "u", "Rjclmm");
+      glvis->Append(u_cf, visspace, "U(psi)", "RjclQmm");
+      glvis->Append(lambda, "lambda", "Rjclmm");
+   }
+   std::unique_ptr<ParaViewDataCollection> paraview_dc;
+   if (paraview)
+   {
+      filename << "r" << ref_levels << "-o" << order;
+      paraview_dc = std::make_unique<ParaViewDataCollection>(filename.str(), &mesh);
+      paraview_dc->SetLevelsOfDetail(order);
+      paraview_dc->SetDataFormat(VTKFormat::BINARY);
+      paraview_dc->SetHighOrderOutput(true);
+      paraview_dc->RegisterField("solution", &x);
+      paraview_dc->SetCycle(0);
+      paraview_dc->SetTime(0.0);
+      paraview_dc->Save();
+   }
+
+   real_t lambda_diff = infinity();
+   for (int i=0; i<100; i++)
+   {
+      alpha = alpha_rule.Get(i);
+      out << "PG iteration " << i + 1 << " with alpha=" << alpha << std::endl;
+      latent_k = latent;
+      latent_k.SetTrueVector();
+
+      solver.Mult(rhs, x_and_latent);
+
+      if (!solver.GetConverged())
+      {
+         out << "Newton Failed to converge in " << solver.GetNumIterations() <<
+             std::endl;
+      }
+      x.SetFromTrueDofs(x_and_latent.GetBlock(0));
+      latent.SetFromTrueDofs(x_and_latent.GetBlock(1));
+
+      if (glvis) { glvis->Update(); }
+      if (paraview_dc)
+      {
+         paraview_dc->SetCycle(i+1);
+         paraview_dc->SetTime(i+1);
+         paraview_dc->Save();
+      }
+
+      subtract(latent, latent_k, lambda);
+      lambda *= 1.0 / pg_functional.GetAlpha();
+
+      if ((lambda_diff = lambda.ComputeL1Error(lambda_prev_cf)) < 1e-8)
+      {
+         out << "  The dual variable, (psi - psi_k)/alpha, converged" << std::endl;
+         out << "PG Converged in " << i + 1
+             << " with final Lambda difference: " << lambda_diff << std::endl;
+         break;
+      }
+      else
+      {
+         out << "  Newton converged in " << solver.GetNumIterations()
+             << " with residual " << solver.GetFinalNorm() << std::endl;
+         out << "  Lambda difference: " << lambda_diff << std::endl;
+      }
+
+      lambda_prev = lambda;
+   }
+   return 0;
+}

miniapps/autodiff/ad_ex6.cpp

@@ -0,0 +1,185 @@
+/// Example 6: Darcy
+#include "mfem.hpp"
+#include "logger.hpp"
+#include "ad_intg.hpp"
+#include "tools.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+struct DarcyFunctional : public ADVectorFunction
+{
+   int dim;
+   // input: q (vector), divq (scalar), u (scalar) -> dim + 1 + 1
+   // output: coefficient for w, divw, v -> dim + 1 + 1 (w, v are test functions)
+   DarcyFunctional(int dim) : ADVectorFunction(dim + 1 + 1, dim + 1 + 1),
+      dim(dim) {}
+   // (q, w) - (div w, u) -> res[w] = q, res[divw] = -u
+   // (div q, v) -> res[v] = div q
+   AD_VEC_IMPL(T, V, M, q_divq_u, res,
+   {
+      res.SetSize(dim + 1 + 1);
+      const V q(q_divq_u.GetData(), dim);
+      const T divq = q_divq_u[dim];
+      const T u = q_divq_u[dim+1];
+
+      V w_cf(res.GetData(), dim);
+      T &divw_cf = res[dim];
+      T &v_cf = res[dim+1];
+      w_cf = q;
+      divw_cf = -u;
+      v_cf = divq;
+   });
+};
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+   MPI_Comm comm = MPI_COMM_WORLD;
+   // file name to be saved
+   std::stringstream filename;
+   filename << "ad-darcy";
+
+   int order = 2;
+   int ref_levels = 3;
+   bool visualization = false;
+   bool paraview = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ref_levels, "-r", "--ref", "Refinement levels");
+   args.AddOption(&visualization, "-vis", "--visualization",
+                  "-no-vis", "--no-visualization",
+                  "Enable visualization, default is false");
+   args.AddOption(&paraview, "-pv", "--paraview",
+                  "-no-pv", "--no-paraview",
+                  "Enable Paraview Export. Default is false");
+   args.ParseCheck();
+   if (myid != 0) { out.Disable(); }
+
+   // Mesh mesh = rhs_fun_circle
+   Mesh ser_mesh = Mesh::MakeCartesian2D(2, 2,
+                                         Element::QUADRILATERAL);
+   const int dim = ser_mesh.Dimension();
+   for (int i = 0; i < ref_levels; i++)
+   {
+      ser_mesh.UniformRefinement();
+   }
+   ParMesh mesh(comm, ser_mesh);
+
+   FunctionCoefficient load_cf([](const Vector &x)
+   {
+      return 2*M_PI * M_PI * std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
+   });
+   DarcyFunctional darcy_functional(dim);
+
+   RT_FECollection flux_fec(order, dim);
+   L2_FECollection potential_fec(order, dim);
+   ParFiniteElementSpace flux_fes(&mesh, &flux_fec);
+   ParFiniteElementSpace potential_fes(&mesh, &potential_fec);
+   QuadratureSpace visspace(&mesh, order+3);
+   const IntegrationRule &ir = IntRules.Get(Geometry::Type::SQUARE, 3*order + 3);
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = flux_fes.GetTrueVSize();
+   offsets[2] = potential_fes.GetTrueVSize();
+   offsets.PartialSum();
+   BlockVector flux_and_potential(offsets);
+   BlockVector rhs(offsets);
+
+   ParGridFunction flux(&flux_fes), potential(&potential_fes);
+
+   flux = 0.0; flux.GetTrueDofs(flux_and_potential.GetBlock(0));
+   potential = 0.0; potential.GetTrueDofs(flux_and_potential.GetBlock(1));
+
+   Array<ParFiniteElementSpace*> fespaces{&flux_fes, &potential_fes};
+   ParBlockNonlinearForm bnlf(fespaces);
+   constexpr ADEval flux_mode = ADEval::VECFE | ADEval::VALUE | ADEval::DIV;
+   constexpr ADEval potential_mode = ADEval::VALUE;
+   bnlf.AddDomainIntegrator(
+      new ADBlockNonlinearFormIntegrator<flux_mode, potential_mode>
+      (darcy_functional, &ir)
+   );
+
+   ParLinearForm b(&potential_fes);
+   b.AddDomainIntegrator(new DomainLFIntegrator(load_cf));
+   b.Assemble();
+   b.ParallelAssemble(rhs.GetBlock(1));
+   rhs.GetBlock(0) = 0.0;
+
+   GMRESSolver lin_solver(comm);
+   lin_solver.SetRelTol(1e-08);
+   lin_solver.SetAbsTol(0.0);
+   lin_solver.SetMaxIter(1e04);
+   lin_solver.SetKDim(100);
+
+   BlockOperator &darcy_op = bnlf.GetGradient(flux_and_potential);
+   Vector Md(flux_fes.GetTrueVSize());
+   HypreParMatrix &M = static_cast<HypreParMatrix&>(darcy_op.GetBlock(0,0));
+   HypreParMatrix &B = static_cast<HypreParMatrix&>(darcy_op.GetBlock(1,0));
+   M.GetDiag(Md);
+   HypreParMatrix invMBt(static_cast<HypreParMatrix&>(darcy_op.GetBlock(0,1)));
+   invMBt.InvScaleRows(Md);
+   std::unique_ptr<HypreParMatrix> S(ParMult(&B, &invMBt));
+   BlockDiagonalPreconditioner prec(offsets);
+   HypreDiagScale invM(M);
+   HypreBoomerAMG invS(*S);
+   invS.SetPrintLevel(0);
+   invM.iterative_mode = false;
+   invS.iterative_mode = false;
+   invS.SetMaxIter(1);
+   prec.SetDiagonalBlock(0, &invM);
+   prec.SetDiagonalBlock(1, &invS);
+   prec.owns_blocks = false;
+
+   lin_solver.SetPreconditioner(prec);
+   lin_solver.SetOperator(darcy_op);
+   lin_solver.Mult(rhs, flux_and_potential);
+   flux.SetFromTrueDofs(flux_and_potential.GetBlock(0));
+   potential.SetFromTrueDofs(flux_and_potential.GetBlock(1));
+
+   if (visualization)
+   {
+      GLVis glvis("localhost", 19916, 400, 350, 3);
+      glvis.Append(flux, "flux", "RjclQmm");
+      glvis.Append(potential, "potential", "Rjclmm");
+   }
+   if (paraview)
+   {
+      std::stringstream pvloc;
+      pvloc << "ParaView/" << filename.str();
+      ParaViewDataCollection paraview_dc(pvloc.str(), &mesh);
+      paraview_dc.SetLevelsOfDetail(order);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.RegisterField("flux", &flux);
+      paraview_dc.RegisterField("potential", &potential);
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetTime(0.0);
+      paraview_dc.Save();
+   }
+
+   FunctionCoefficient exact_potential([](const Vector &x)
+   {
+      return std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
+   });
+   VectorFunctionCoefficient exact_flux(dim, [](const Vector &x, Vector &q)
+   {
+      // flux = - grad u
+      q.SetSize(x.Size());
+      q[0] = -M_PI*std::cos(M_PI*x[0])*std::sin(M_PI*x[1]);
+      q[1] = -M_PI*std::sin(M_PI*x[0])*std::cos(M_PI*x[1]);
+   });
+   out << "L2 Error in Potential: "
+       << potential.ComputeL2Error(exact_potential) << std::endl;
+   out << "L2 Error in Flux: "
+       << flux.ComputeL2Error(exact_flux) << std::endl;
+
+   return 0;
+}

miniapps/autodiff/ad_intg.hpp

@@ -0,0 +1,334 @@
+/// Templated AD (block) nonlinear form integrators definitions
+#pragma once
+#include "mfem.hpp"
+#include "ad_native.hpp"
+
+namespace mfem
+{
+// A helper struct to pass the loop index as a template parameter
+template<std::size_t N>
+struct __loop_index
+{
+   static const constexpr size_t value = N;
+   constexpr operator std::size_t() const { return N; }
+};
+// loop over indeces at compile time
+template <class F, std::size_t... Is>
+void _constexpr_for(F func, std::index_sequence<Is...>)
+{
+   (func(__loop_index<Is> {}), ...);
+}
+// @brief ADQuadEvalMode is an enumeration for the evaluation modes of the ADEnergy class.
+// For example, if you want to evaluate the value and gradient of the function, you can use
+// constexpr auto mode = ADEval::VALUE | ADEval::GRAD;
+enum class ADEval
+{
+   QVALUE  = 1 << 0, // u(T, ip) (quadrature value)
+   VALUE   = 1 << 1, // u(T, ip)
+   GRAD    = 1 << 2, // grad u(T, ip)
+   DIV     = 1 << 3, // div u(T, ip) (not yet implemented)
+   CURL    = 1 << 4, // curl u(T, ip) (not yet implemented)
+   Hessian = 1 << 5, // D^2 u(T, ip) (not yet implemented)
+
+   VECTOR  = 1 << 6, // vector-valued scalar FE
+   VECFE   = 1 << 7, // vector-valued vector FE (not yet implemented)
+   NUMOPT  = 1 << 8, // number of options. If change options, change this value to last
+};
+
+constexpr ADEval operator|(ADEval a, ADEval b)
+{
+   return static_cast<ADEval>(static_cast<int>(a) | static_cast<int>(b));
+}
+constexpr ADEval operator&(ADEval a, ADEval b)
+{
+   return static_cast<ADEval>(static_cast<int>(a) & static_cast<int>(b));
+}
+inline constexpr ADEval operator~(ADEval mode)
+{
+   return static_cast<ADEval>(~static_cast<int>(mode));
+}
+inline constexpr bool hasFlag(ADEval mode, ADEval flag)
+{
+   return (mode & flag) == flag;
+}
+
+template <ADEval mode>
+constexpr bool isValidADEval()
+{
+   constexpr auto INVALID = ADEval::Hessian;
+   if constexpr (static_cast<int>(mode & INVALID) != 0) { return false; }
+   if constexpr (hasFlag(mode, ADEval::QVALUE))
+   {
+      // QVALUE cannot be combined with other modes except VECTOR
+      return static_cast<int>(mode & (~(ADEval::QVALUE | ADEval::VECTOR))) == 0;
+   }
+   if constexpr (hasFlag(mode, ADEval::VECFE))
+   {
+      return !hasFlag(mode,
+                      ADEval::VECTOR); // VECTOR is only for vector-valued scalar FE
+   }
+   return true;
+}
+
+template <ADEval... modes>
+class ADBlockNonlinearFormIntegrator;
+
+template <ADEval mode>
+class ADNonlinearFormIntegrator : public NonlinearFormIntegrator
+{
+   static_assert(isValidADEval<mode>(),
+                 "ADNonlinearFormIntegrator: Invalid ADEval mode");
+protected:
+   ADFunction &f;
+
+private:
+   Vector x, jac;
+   DenseMatrix H, Hx;
+
+   // only if ADEvalInput::VECTOR. Each column corresponds to a vector component
+   DenseMatrix xmat, jacMat, Hs, Hxsub;
+   DenseMatrix elfun_matview, elvectmat, partelmat;
+
+   DenseMatrix allshapes; // all shapes, [?shape, ?dshape]
+   Vector shape, shape1, shape2;
+   DenseMatrix vshape, vshape1, vshape2;
+   DenseMatrix gshape, gshape1, gshape2;
+   Vector divshape, divshape1, divshape2;
+   DenseMatrix curlshape, curlshape1, curlshape2;
+   Vector nor;
+   // DenseMatrix d2shape, d2shape1, d2shape2; // for hessian. Not implemented yet.
+public:
+   ADNonlinearFormIntegrator(ADFunction &f, IntegrationRule *ir = nullptr)
+      : NonlinearFormIntegrator(ir), f(f) {}
+
+   const IntegrationRule* GetDefaultIntegrationRule(
+      const FiniteElement& trial_fe, const FiniteElement& test_fe,
+      const ElementTransformation& trans) const override
+   {
+      int order = std::max(trial_fe.GetOrder(), test_fe.GetOrder());
+      return &IntRules.Get(trans.GetGeometryType(), order*2 + 2);
+   }
+
+   /// Compute the local energy
+   real_t GetElementEnergy(const FiniteElement &el,
+                           ElementTransformation &Tr,
+                           const Vector &elfun) override;
+
+   /// Perform the local action of the NonlinearFormIntegrator
+   void AssembleElementVector(const FiniteElement &el,
+                              ElementTransformation &Tr,
+                              const Vector &elfun, Vector &elvect) override;
+
+   /// Assemble the local gradient matrix
+   void AssembleElementGrad(const FiniteElement &el,
+                            ElementTransformation &Tr,
+                            const Vector &elfun, DenseMatrix &elmat) override;
+
+   /// @brief Perform the local action of the NonlinearFormIntegrator resulting
+   /// from a face integral term.
+   void AssembleFaceVector(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Tr,
+                           const Vector &elfun, Vector &elvect) override;
+
+
+   /// @brief Assemble the local action of the gradient of the
+   /// NonlinearFormIntegrator resulting from a face integral term.
+   void AssembleFaceGrad(const FiniteElement &el1,
+                         const FiniteElement &el2,
+                         FaceElementTransformations &Tr,
+                         const Vector &elfun, DenseMatrix &elmat) override;
+
+
+protected:
+
+   // Initialize shapes to [?value_shapes, ?grad_shapes]
+   // and make value_shapes and grad_shapes reference to
+   // allshapes.
+   inline int InitInputShapes(const FiniteElement &el,
+                              ElementTransformation &Tr,
+                              DenseMatrix &shapes);
+
+   // Calculate parameter, shape, dshape at the given integration point
+   inline void CalcInputShapes(const FiniteElement &el,
+                               ElementTransformation &Tr,
+                               const IntegrationPoint &ip,
+                               DenseMatrix &shapes);
+   template <ADEval... modes>
+   friend class ADBlockNonlinearFormIntegrator;
+private:
+};
+
+template <ADEval... modes>
+class ADBlockNonlinearFormIntegrator : public BlockNonlinearFormIntegrator
+{
+public:
+   const IntegrationRule *IntRule = nullptr;
+
+protected:
+   constexpr static size_t numSpaces = sizeof...(modes);
+   static constexpr std::array<ADEval, sizeof...(modes)> modes_arr = {modes...};
+   ADFunction &f;
+   const IntegrationRule* GetIntegrationRule(
+      const FiniteElement& trial_fe, const FiniteElement& test_fe,
+      const ElementTransformation& trans) const;
+
+   /** @brief Returns an integration rule based on the arguments and
+              internal state. (Version for identical trial_fe and test_fe)
+
+       @see GetIntegrationRule(const FiniteElement*, const FiniteElement*,
+            const ElementTransformation*)
+   */
+   const IntegrationRule* GetIntegrationRule(
+      const FiniteElement& el,
+      const ElementTransformation& trans) const;
+
+private:
+   Array<int> vdim;
+   Vector x, jac;
+   std::vector<Vector> xvar, jacVar;
+   DenseMatrix H;
+   DenseMatrix Hsub;
+   DenseMatrix Hx;
+   DenseMatrix Hxsub;
+
+   // only if ADEvalInput::VECTOR. Each column corresponds to a vector component
+   std::vector<DenseMatrix> xmat, jacVarMat, Hs;
+   std::vector<DenseMatrix> elfun_matview, elvectmat, partelmat;
+
+   std::vector<DenseMatrix> allshapes; // all shapes, [?shape, ?dshape]
+   std::vector<Vector> shape, shape1, shape2;
+   std::vector<DenseMatrix> vshape, vshape1, vshape2;
+   std::vector<DenseMatrix> gshape, gshape1, gshape2;
+   std::vector<Vector> divshape, divshape1, divshape2;
+   std::vector<DenseMatrix> curlshape, curlgshape1, curlgshape2;
+   Vector nor;
+   // DenseMatrix d2shape, d2shape1, d2shape2; // for hessian. Not implemented yet.
+public:
+   ADBlockNonlinearFormIntegrator(ADFunction &f,
+                                  const IntegrationRule *ir = nullptr)
+      : IntRule(ir), f(f), vdim(numSpaces)
+      , allshapes(numSpaces)
+      , xvar(numSpaces), jacVar(numSpaces)
+      , Hx(numSpaces)
+      , xmat(numSpaces), jacVarMat(numSpaces)
+      , Hs(numSpaces), Hxsub(numSpaces)
+      , elfun_matview(numSpaces), elvectmat(numSpaces)
+      , partelmat(numSpaces)
+      , shape(numSpaces), shape1(numSpaces), shape2(numSpaces)
+      , vshape(numSpaces), vshape1(numSpaces), vshape2(numSpaces)
+      , gshape(numSpaces), gshape1(numSpaces), gshape2(numSpaces)
+      , divshape(numSpaces), divshape1(numSpaces), divshape2(numSpaces)
+      , curlshape(numSpaces), curlgshape1(numSpaces), curlgshape2(numSpaces)
+   { vdim = 1; }
+
+   ADBlockNonlinearFormIntegrator(ADFunction &f, std::initializer_list<int> vdim,
+                                  const IntegrationRule *ir = nullptr)
+      : ADBlockNonlinearFormIntegrator(f, ir), vdim(vdim)
+   {}
+
+   virtual void SetIntRule(const IntegrationRule *ir)
+   { IntRule = ir; }
+
+   /** @brief Prescribe a fixed IntegrationRule to use. Sets the NURBS patch
+              integration rule to null.
+
+       @see SetIntRule(const IntegrationRule*)
+   */
+   void SetIntegrationRule(const IntegrationRule &ir) { SetIntRule(&ir); }
+
+   /** @brief Directly return the IntRule pointer (possibly null) without
+       checking for NURBS patch rules or falling back on a default. */
+   const IntegrationRule *GetIntRule() const { return IntRule; }
+
+   /** @brief Equivalent to GetIntRule, but retained for backward
+       compatibility with applications. */
+   const IntegrationRule *GetIntegrationRule() const { return GetIntRule(); }
+
+
+   /// Compute the local energy
+   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
+                           ElementTransformation &Tr,
+                           const Array<const Vector*> &elfun) override;
+
+   /// Perform the local action of the NonlinearFormIntegrator
+   void AssembleElementVector(const Array<const FiniteElement *>&el,
+                              ElementTransformation &Tr,
+                              const Array<const Vector *>&elfun,
+                              const Array<Vector *>&elvect) override;
+
+   /// Assemble the local gradient matrix
+   void AssembleElementGrad(const Array<const FiniteElement *>&el,
+                            ElementTransformation &Tr,
+                            const Array<const Vector *>&elfun,
+                            const Array2D<DenseMatrix *>&elmat) override;
+
+   /// @brief Perform the local action of the NonlinearFormIntegrator resulting
+   /// from a face integral term.
+   void AssembleFaceVector(const Array<const FiniteElement *>&el1,
+                           const Array<const FiniteElement *>&el2,
+                           FaceElementTransformations &Tr,
+                           const Array<const Vector *>&elfun,
+                           const Array<Vector *>&elvect) override;
+
+
+   /// @brief Assemble the local action of the gradient of the
+   /// NonlinearFormIntegrator resulting from a face integral term.
+   void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
+                         const Array<const FiniteElement *>&el2,
+                         FaceElementTransformations &Tr,
+                         const Array<const Vector *>&elfun,
+                         const Array2D<DenseMatrix *>&elmat) override;
+
+
+protected:
+
+   const IntegrationRule* GetIntegrationRule(
+      const Array<const FiniteElement *>& trial_fe,
+      const Array<const FiniteElement *>& test_fe,
+      const ElementTransformation& trans) const
+   {
+      if (IntRule) { return IntRule; }
+      return GetDefaultIntegrationRule(trial_fe, test_fe, trans);
+   }
+
+   const IntegrationRule* GetIntegrationRule(
+      const Array<const FiniteElement *>& el,
+      const ElementTransformation& trans) const
+   {
+      if (IntRule) { return IntRule; }
+      return GetDefaultIntegrationRule(el, el, trans);
+   }
+
+   virtual const IntegrationRule* GetDefaultIntegrationRule(
+      const Array<const FiniteElement *>& trial_fe,
+      const Array<const FiniteElement *>& test_fe,
+      const ElementTransformation& trans) const
+   {
+      int order = 0;
+      for (int i=0; i<trial_fe.Size(); i++)
+      {
+         order = std::max(order, trial_fe[i]->GetOrder());
+      }
+      for (int i=0; i<test_fe.Size(); i++)
+      {
+         order = std::max(order, test_fe[i]->GetOrder());
+      }
+      return &IntRules.Get(trans.GetGeometryType(), order*2 + 2);
+   }
+
+   std::array<int, sizeof...(modes)> InitInputShapes(
+      const Array<const FiniteElement *>& el,
+      ElementTransformation &Tr,
+      std::vector<DenseMatrix> &shapes);
+
+   void CalcInputShapes(
+      const Array<const FiniteElement *>& el,
+      ElementTransformation &Tr,
+      const IntegrationPoint &ip,
+      std::vector<DenseMatrix> &allshapes);
+
+private:
+};
+} // namespace mfem
+#include "_ad_intg.hpp"

miniapps/autodiff/ad_native.cpp

@@ -0,0 +1,327 @@
+#include "ad_native.hpp"
+
+namespace mfem
+{
+int Evaluator::GetSize(const param_t &param)
+{
+   return std::visit([](auto arg)
+   {
+      using T = std::decay_t<decltype(arg)>;
+      if constexpr (std::is_same_v<T, real_t>)
+      {
+         return 1;
+      }
+      if constexpr (std::is_same_v<T, Vector>)
+      {
+         return arg.Size();
+      }
+      if constexpr (std::is_same_v<T, DenseMatrix>)
+      {
+         return arg.TotalSize();
+      }
+      if constexpr (std::is_same_v<T, const real_t*>)
+      {
+         return 1;
+      }
+      if constexpr (std::is_same_v<T, const Vector*>)
+      {
+         return arg->Size();
+      }
+      if constexpr (std::is_same_v<T, const DenseMatrix*>)
+      {
+         return arg->Height()*arg->Width();
+      }
+      if constexpr (std::is_same_v<T, Coefficient*>)
+      {
+         return 1;
+      }
+      if constexpr (std::is_same_v<T, VectorCoefficient*>)
+      {
+         return arg->GetVDim();
+      }
+      if constexpr (std::is_same_v<T, MatrixCoefficient*>)
+      {
+         return arg->GetHeight() * arg->GetWidth();
+      }
+      if constexpr (std::is_same_v<T, const GridFunction*>)
+      {
+         return arg->FESpace()->GetVDim();
+      }
+      if constexpr (std::is_same_v<T, const QuadratureFunction*>)
+      {
+         return arg->GetVDim();
+      }
+      MFEM_ABORT("Evaluator: Unsupported parameter type");
+      return 0;
+   }, param);
+}
+
+Evaluator::~Evaluator()
+{
+   for (int i=0; i<params.size(); i++)
+   {
+      if (owns[i])
+      {
+         std::visit([](auto &arg)
+         {
+            using T = std::decay_t<decltype(arg)>;
+            if constexpr (std::is_pointer_v<T>)
+            {
+               delete arg;
+            }
+         }, params[i]);
+      }
+   }
+}
+int Evaluator::Add(param_t param, bool eval_owns)
+{
+   int idx = params.size();
+   params.push_back(param);
+   offsets.Append(offsets.Last() + GetSize(param));
+   val.Update(offsets);
+   owns.Append(eval_owns);
+   std::visit([&](auto arg)
+   {
+      using T = std::decay_t<decltype(arg)>;
+      if constexpr (std::is_same_v<T, real_t>)
+      {
+         MFEM_VERIFY(eval_owns==false,
+                     "Evaluator::Add: real_t parameter cannot own the value");
+         val.GetBlock(idx) = arg;
+      }
+      if constexpr (std::is_same_v<T, Vector>)
+      {
+         MFEM_VERIFY(eval_owns==false,
+                     "Evaluator::Add: real_t parameter cannot own the value");
+         val.GetBlock(idx) = arg;
+      }
+      if constexpr (std::is_same_v<T, DenseMatrix>)
+      {
+         MFEM_VERIFY(eval_owns==false,
+                     "Evaluator::Add: real_t parameter cannot own the value");
+         Vector v(arg.GetData(), arg.TotalSize());
+         val.GetBlock(idx) = v;
+      }
+   }, param);
+   return idx;
+}
+void Evaluator::Replace(size_t i, param_t param)
+{
+   MFEM_VERIFY(i < params.size(),
+               "Evaluator::Set: index out of range");
+   params[i] = param;
+   int size = GetSize(param);
+   MFEM_VERIFY(size == offsets[i+1] - offsets[i],
+               "Evaluator::Set: size mismatch for parameter at index " << i
+               << ": expected " << (offsets[i+1] - offsets[i]) << ", got " << size);
+}
+
+const Vector& Evaluator::Eval(int i, ElementTransformation &Tr,
+                              const IntegrationPoint &ip) const
+{
+   std::visit([&](auto arg)
+   {
+      Vector &v = this->val.GetBlock(i);
+      using T = std::decay_t<decltype(arg)>;
+      if constexpr (std::is_same_v<T, real_t> ||
+                    std::is_same_v<T, Vector> ||
+                    std::is_same_v<T, DenseMatrix>)
+      {
+         // Already stored, do nothing
+         return;
+      }
+      if constexpr (std::is_same_v<T, const real_t*>)
+      {
+         v = *arg;
+         return;
+      }
+      if constexpr (std::is_same_v<T, const Vector*>)
+      {
+         v = *arg;
+         return;
+      }
+      if constexpr (std::is_same_v<T, const DenseMatrix*>)
+      {
+         DenseMatrix m(v.GetData(), arg->Height(), arg->Width());
+         m = *arg;
+         return;
+      }
+      if constexpr (std::is_same_v<T, Coefficient*>)
+      {
+         v(0) = arg->Eval(Tr, ip);
+         return;
+      }
+      if constexpr (std::is_same_v<T, VectorCoefficient*>)
+      {
+         arg->Eval(v, Tr, ip);
+         return;
+      }
+      if constexpr (std::is_same_v<T, MatrixCoefficient*>)
+      {
+         DenseMatrix m(v.GetData(), arg->GetHeight(), arg->GetWidth());
+         arg->Eval(m, Tr, ip);
+         return;
+      }
+      if constexpr (std::is_same_v<T, const GridFunction*>)
+      {
+         arg->GetVectorValue(Tr, ip, v);
+         return;
+      }
+      if constexpr (std::is_same_v<T, const QuadratureFunction*>)
+      {
+         arg->GetValues(Tr.ElementNo, ip.index, v);
+         return;
+      }
+   }, params[i]);
+   return this->val.GetBlock(i);
+}
+
+void ADFunction::Gradient(const Vector &x, ElementTransformation &Tr,
+                          const IntegrationPoint &ip,
+                          Vector &J) const
+{
+   ProcessParameters(Tr, ip);
+   Gradient(x, J);
+}
+void ADFunction::Gradient(const Vector &x, Vector &J) const
+{
+   MFEM_ASSERT(x.Size() == n_input,
+               "ADFunction::Gradient: x.Size() must match n_input");
+   J.SetSize(x.Size());
+   ADVector x_ad(x);
+   for (int i=0; i < n_input; i++)
+   {
+      x_ad[i].gradient = 1.0;
+      ADReal_t result = (*this)(x_ad);
+      J[i] = result.gradient;
+      x_ad[i].gradient = 0.0;
+   }
+}
+
+void ADFunction::Hessian(const Vector &x, ElementTransformation &Tr,
+                         const IntegrationPoint &ip,
+                         DenseMatrix &H) const
+{
+   ProcessParameters(Tr, ip);
+   Hessian(x, H);
+}
+
+void ADFunction::Hessian(const Vector &x, DenseMatrix &H) const
+{
+   MFEM_ASSERT(x.Size() == n_input,
+               "ADFunction::Hessian: x.Size() must match n_input");
+   H.SetSize(x.Size(), x.Size());
+   AD2Vector x_ad(x);
+   for (int i=0; i<n_input; i++) // Loop for the first derivative
+   {
+      x_ad[i].value.gradient = 1.0;
+      for (int j=0; j<=i; j++)
+      {
+         x_ad[j].gradient.value = 1.0;
+         AD2Real_t result = (*this)(x_ad);
+         H(j, i) = result.gradient.gradient;
+         H(i, j) = result.gradient.gradient;
+         x_ad[j].gradient.value = 0.0; // Reset gradient for next iteration
+      }
+      x_ad[i].value.gradient = 0.0;
+   }
+}
+
+void ADVectorFunction::Gradient(const Vector &x, DenseMatrix &J) const
+{
+   MFEM_ASSERT(x.Size() == n_input,
+               "ADVectorFunction::Gradient: x.Size() must match n_input");
+   ADVector x_ad(x);
+   ADVector Fx(n_output);
+   J.SetSize(n_output, n_input);
+   for (int i=0; i<n_input; i++)
+   {
+      x_ad[i].gradient = 1.0;
+      Fx = ADReal_t();
+      (*this)(x_ad, Fx);
+      for (int j=0; j<n_output; j++)
+      {
+         J(j,i) = Fx[j].gradient;
+      }
+      x_ad[i].gradient = 0.0; // Reset gradient for next iteration
+   }
+}
+
+void ADVectorFunction::Hessian(const Vector &x, DenseTensor &H) const
+{
+   MFEM_ASSERT(x.Size() == n_input,
+               "ADVectorFunction::Gradient: x.Size() must match n_input");
+   AD2Vector x_ad(x);
+   AD2Vector Fx(n_output);
+   H.SetSize(n_input, n_input, n_output);
+   for (int i=0; i<n_input; i++) // Loop for the first derivative
+   {
+      x_ad[i].value.gradient = 1.0;
+      for (int j=0; j<=i; j++)
+      {
+         x_ad[j].gradient.value = 1.0;
+         Fx = AD2Real_t();
+         (*this)(x_ad, Fx);
+         for (int k=0; k<n_output; k++)
+         {
+            H(j, i, k) = Fx[k].gradient.gradient;
+            H(i, j, k) = Fx[k].gradient.gradient;
+         }
+         x_ad[j].gradient.value = 0.0; // Reset gradient for next iteration
+      }
+      x_ad[i].value.gradient = 0.0;
+   }
+}
+
+
+Lagrangian Lagrangian::AddEqConstraint(ADFunction &constraint,
+                                       real_t target)
+{
+   eq_con.push_back(&constraint);
+   int numCon = eq_con.size();
+   eq_rhs.SetSize(numCon);
+   eq_rhs[numCon - 1] = target;
+   return *this;
+}
+
+void Lagrangian::ProcessParameters(ElementTransformation &Tr,
+                                   const IntegrationPoint &ip) const
+{
+   objective.ProcessParameters(Tr, ip);
+   for (auto *con : eq_con) { con->ProcessParameters(Tr, ip); }
+}
+
+ALFunctional ALFunctional::AddEqConstraint(ADFunction &constraint,
+                                           real_t target)
+{
+   eq_con.push_back(&constraint);
+   int numCon = eq_con.size();
+   eq_rhs.SetSize(numCon);
+   lambda.SetSize(numCon);
+
+   eq_rhs[numCon - 1] = target;
+   lambda[numCon - 1] = 0.0;
+   return *this;
+}
+
+void ALFunctional::SetLambda(const Vector &lambda)
+{
+   MFEM_VERIFY(lambda.Size() == this->lambda.Size(),
+               "ALFunctional: lambda size mismatch");
+   this->lambda = lambda;
+}
+
+void ALFunctional::SetPenalty(real_t mu)
+{
+   MFEM_VERIFY(mu >= 0.0, "ALFunctional: mu must be non-negative");
+   this->penalty = mu;
+}
+
+void ALFunctional::ProcessParameters(ElementTransformation &Tr,
+                                     const IntegrationPoint &ip) const
+{
+   objective.ProcessParameters(Tr, ip);
+   for (auto *con : eq_con) { con->ProcessParameters(Tr, ip); }
+}
+
+} // namespace mfem

miniapps/autodiff/ad_native.hpp

@@ -0,0 +1,829 @@
+#pragma once
+
+#include "mfem.hpp"
+#include "tadvector.hpp"
+#include "taddensemat.hpp"
+
+namespace mfem
+{
+template <typename T>
+constexpr auto type_name()
+{
+#if defined(__clang__)
+   return std::string_view(__PRETTY_FUNCTION__);
+#elif defined(__GNUC__)
+   return std::string_view(__PRETTY_FUNCTION__);
+#elif defined(_MSC_VER)
+   return std::string_view(__FUNCSIG__);
+#else
+   return std::string_view("unknown");
+#endif
+}
+template <typename value_type, typename gradient_type, typename other_type>
+MFEM_HOST_DEVICE
+inline future::dual<value_type, gradient_type> max(
+   future::dual<value_type, gradient_type> a, other_type b);
+
+inline real_t max(const real_t a, const real_t b) { return std::max(a,b); }
+
+template <typename value_type, typename gradient_type, typename other_type>
+MFEM_HOST_DEVICE
+inline future::dual<value_type, gradient_type> min(
+   future::dual<value_type, gradient_type> a, other_type b);
+
+MFEM_HOST_DEVICE
+inline real_t min(const real_t a, const real_t b) { return std::min(a,b); }
+
+// Use mfem-native autodiff types
+// If other autodiff libraries are used,
+// define ADReal_t, ADVector, ADMatrix, ... types accordingly.
+
+// First order dual
+typedef future::dual<real_t, real_t> ADReal_t;
+typedef TAutoDiffVector<ADReal_t> ADVector;
+typedef TAutoDiffDenseMatrix<ADReal_t> ADMatrix;
+
+// second order dual (nested dual)
+typedef future::dual<ADReal_t, ADReal_t> AD2Real_t;
+typedef TAutoDiffVector<AD2Real_t> AD2Vector;
+typedef TAutoDiffDenseMatrix<AD2Real_t> AD2Matrix;
+
+class Evaluator
+{
+   // To add a new parameter type,
+   // implement GetSize() and Eval() method
+public:
+   using param_t = std::variant<
+                   real_t, Vector, DenseMatrix, // pass by value
+                   const real_t*, const Vector*, const DenseMatrix*, // pass by pointer
+                   Coefficient*, VectorCoefficient*, MatrixCoefficient*,
+                   const GridFunction*,
+                   const QuadratureFunction*>;
+private:
+   Array<int> offsets;
+   std::vector<param_t> params;
+
+   mutable Vector loc_vec_val;
+   mutable DenseMatrix loc_mat_val;
+
+public:
+   mutable BlockVector val;
+   mutable Array<bool> owns;
+   Evaluator(): offsets{0} {}
+   Evaluator(int capacity)
+      : offsets{0}
+   {
+      val.SetSize(capacity);
+      val.SetSize(0);
+   }
+   virtual ~Evaluator();
+   // Add a parameter to the evaluator
+   int Add(param_t param, bool eval_owns = false);
+   int Add(Vector &v)
+   {
+      if (dynamic_cast<GridFunction*>(&v))
+      {
+         MFEM_WARNING("Adding GridFunction by value, instead of its pointer. "
+                      "This result in the whole GridFunction value will be used at each quadrature point, "
+                      "which is likely not what you want. "
+                      "Use Add(const GridFunction*) instead.");
+      }
+      return Add((param_t)v);
+   }
+   int Add(const Vector &v)
+   {
+      if (dynamic_cast<const GridFunction*>(&v))
+      {
+         MFEM_WARNING("Adding GridFunction by value, instead of its pointer. "
+                      "This result in the whole GridFunction value will be used at each quadrature point, "
+                      "which is likely not what you want. "
+                      "Use Add(const GridFunction*) instead.");
+      }
+      return Add((param_t)v);
+   }
+   // Replace a parameter at index i with a new parameter
+   // The output size of param should match the size of the old parameter
+   void Replace(size_t i, param_t param);
+   param_t Get(size_t i) const
+   {
+      MFEM_VERIFY(i >= 0 && i < params.size(),
+                  "Evaluator::Get: index out of range");
+      return params[i];
+   }
+
+   // Evaluate all parameters at once
+   // and return the block vector
+   const BlockVector &Eval(ElementTransformation &Tr,
+                           const IntegrationPoint &ip) const
+   {
+      for (int i=0; i<params.size(); i++)
+      { Eval(i, Tr, ip); }
+      return val;
+   }
+
+   // Evaluate the parameter at index i
+   // this will update the val block vector, and return the corresponding block
+   const Vector& Eval(int i, ElementTransformation &Tr,
+                      const IntegrationPoint &ip) const;
+   static int GetSize(const param_t &param);
+   int GetSize(size_t i) const
+   {
+      MFEM_VERIFY(i >= 0 && i < offsets.Size() - 1,
+                  "Evaluator::GetSize: index out of range");
+      return offsets[i+1] - offsets[i];
+   }
+   void Project(QuadratureFunction &qf)
+   {
+      const int vdim = offsets.Last();
+      qf.SetVDim(vdim);
+
+      QuadratureSpaceBase &qspace = *qf.GetSpace();
+      Vector qf_view(qf.GetData(), vdim);
+      for (int i=0; i<qspace.GetNE(); i++)
+      {
+         ElementTransformation &Tr = *qspace.GetTransformation(i);
+         const IntegrationRule &ir = qspace.GetIntRule(i);
+         for (int j=0; j<ir.GetNPoints(); j++)
+         {
+            const IntegrationPoint &ip = ir.IntPoint(j);
+            qf_view = Eval(Tr, ip);
+            qf_view.SetData(qf.GetData() + vdim);
+         }
+      }
+      return;
+   }
+   int GetVDim() const { return offsets.Last(); }
+};
+
+class EvaluatorCF : public Coefficient
+{
+   Evaluator &evaluator;
+   int idx;
+   const real_t &val;
+public:
+   EvaluatorCF(Evaluator &evaluator_, int outer_idx=0, int inner_idx=0)
+      : evaluator(evaluator_)
+      , idx(outer_idx)
+      , val(evaluator.val.GetBlock(outer_idx)(inner_idx)) {}
+   real_t Eval(ElementTransformation &Tr, const IntegrationPoint &ip) override
+   {
+      evaluator.Eval(idx, Tr, ip);
+      return val;
+   }
+};
+class EvaluatorVCF : public VectorCoefficient
+{
+   Evaluator &evaluator;
+   int idx;
+public:
+   EvaluatorVCF(Evaluator &evaluator, int idx=-1)
+      : VectorCoefficient(idx == -1 ? evaluator.GetVDim() :
+                          evaluator.val.GetBlock(idx).Size())
+      , evaluator(evaluator)
+      , idx(idx)
+   { }
+
+   void Eval(Vector &V, ElementTransformation &Tr,
+             const IntegrationPoint &ip) override
+   {
+      if (idx == -1) { V = evaluator.Eval(Tr, ip); }
+      else { V = evaluator.Eval(idx, Tr, ip); }
+   }
+};
+class EvaluatorMCF : public MatrixCoefficient
+{
+   Evaluator &evaluator;
+   int idx;
+   const DenseMatrix val;
+public:
+   EvaluatorMCF(Evaluator &evaluator, int h, int w, int idx=0)
+      : MatrixCoefficient(h, w)
+      , evaluator(evaluator)
+      , idx(idx)
+      , val(evaluator.val.GetBlock(idx).GetData(), h, w)
+   {
+      MFEM_VERIFY(evaluator.val.GetBlock(idx).Size() == h*w,
+                  "EvaluatorMCF: size mismatch");
+   }
+   void Eval(DenseMatrix &M, ElementTransformation &Tr,
+             const IntegrationPoint &ip) override
+   {
+      evaluator.Eval(idx, Tr, ip);
+      M = val;
+   }
+};
+
+class ADFunction
+{
+protected:
+
+   int AddParameter(Evaluator::param_t param)
+   { return evaluator.Add(param); }
+
+   void ReplaceParameter(int i, Evaluator::param_t param)
+   { evaluator.Replace(i, param); }
+
+   Evaluator evaluator;
+public:
+   virtual void ProcessParameters(ElementTransformation &Tr,
+                                  const IntegrationPoint &ip) const
+   { ProcessParameters(evaluator.Eval(Tr, ip)); }
+   virtual void ProcessParameters(const BlockVector &param_val) const
+   { }
+
+   const int n_input;
+   ADFunction(int n_input): n_input(n_input) {}
+   // Constructor with capacity for evaluator.
+   // This is useful when the parameter size is known in advance,
+   // so that we can get references to the parameters at construction time.
+   ADFunction(int n_input, int capacity)
+      : n_input(n_input), evaluator(capacity)
+   {
+      MFEM_ASSERT(n_input > 0, "ADFunction: n_input must be positive");
+   }
+   // default evaluator
+   virtual real_t operator()(const Vector &x) const
+   { MFEM_ABORT("Not implemented. Use AD_IMPL macro to implement all path"); }
+   virtual real_t operator()(const Vector &x, ElementTransformation &Tr,
+                             const IntegrationPoint &ip) const
+   { ProcessParameters(Tr, ip); return (*this)(x); }
+
+   // default Jacobian evaluator
+   virtual ADReal_t operator()(const ADVector &x) const
+   { MFEM_ABORT("Not implemented. Use MAKE_AD_FUNCTOR macro to create derived structure"); }
+
+   // default Hessian evaluator
+   virtual AD2Real_t operator()(const AD2Vector &x) const
+   { MFEM_ABORT("Not implemented. Use MAKE_AD_FUNCTOR macro to create derived structure"); }
+
+   // Evaluate the gradient, using forward mode autodiff
+   virtual void Gradient(const Vector &x, ElementTransformation &Tr,
+                         const IntegrationPoint &ip, Vector &J) const;
+   virtual void Gradient(const Vector &x, Vector &J) const;
+   // Evaluate the Hessian, using forward over forward autodiff
+   // The Hessian assumed to be symmetric.
+   virtual void Hessian(const Vector &x, ElementTransformation &Tr,
+                        const IntegrationPoint &ip,
+                        DenseMatrix &H) const;
+   virtual void Hessian(const Vector &x, DenseMatrix &H) const;
+};
+
+// We currently only support Jacobian.
+// To consistent with ADFunction, which returns
+// evaluate: scalar, Gradient: vector, Hessian: matrix,
+// we overrode the Gradient for evaulation, and Hessian for Jacobian
+// To be used with ADNonlinearFormIntegrator or ADBlockNonlinearFormIntegrator,
+// n_input and n_output must be the same.
+struct ADVectorFunction : public ADFunction
+{
+
+   int n_output;
+   ADVectorFunction(int n_input, int n_output)
+      : ADFunction(n_input), n_output(n_output)
+   {
+      MFEM_ASSERT(n_input > 0 && n_output > 0,
+                  "ADVectorFunction: n_input and n_output must be positive");
+   }
+
+   void operator()(const Vector &x, ElementTransformation &Tr,
+                   const IntegrationPoint &ip,
+                   Vector &F) const
+   { ProcessParameters(Tr, ip); (*this)(x, F); }
+
+   // Derived struct should implement the following methods.
+   // Use AD_VEC_IMPL macro to implement them.
+   virtual void operator()(const Vector &x, Vector &F) const = 0;
+   virtual void operator()(const ADVector &x, ADVector &F) const = 0;
+   virtual void operator()(const AD2Vector &x, AD2Vector &F) const = 0;
+
+   void Gradient(const Vector &x, ElementTransformation &Tr,
+                 const IntegrationPoint &ip, DenseMatrix &J) const
+   { ProcessParameters(Tr, ip); Gradient(x, J); }
+
+   void Gradient(const Vector &x, DenseMatrix &J) const;
+
+   void Hessian(const Vector &x, ElementTransformation &Tr,
+                const IntegrationPoint &ip,
+                DenseTensor &H) const
+   { ProcessParameters(Tr, ip); Hessian(x, H); }
+
+   void Hessian(const Vector &x, DenseTensor &H) const;
+
+   // To support ADNonlinearFormIntegrator and ADVectorNonlinearFormIntegrator
+   void Gradient(const Vector &x, ElementTransformation &Tr,
+                 const IntegrationPoint &ip, Vector &F) const override final
+   { (*this)(x, Tr, ip, F); }
+
+   void Gradient(const Vector &x, Vector &F) const override final
+   { (*this)(x, F); }
+
+   // To support ADNonlinearFormIntegrator and ADVectorNonlinearFormIntegrator
+   void Hessian(const Vector &x, ElementTransformation &Tr,
+                const IntegrationPoint &ip,
+                DenseMatrix &J) const override final
+   { this->Gradient(x, Tr, ip, J); }
+
+   void Hessian(const Vector &x, DenseMatrix &J) const override final
+   { this->Gradient(x, J); }
+
+   real_t operator()(const Vector &x) const override final
+   {
+      MFEM_ABORT("ADVectorFunction::operator(): This method should not be called. "
+                 "Use ADVectorFunction::operator(const Vector &x, Vector &F) instead.");
+   }
+   ADReal_t operator()(const ADVector &x) const override final
+   {
+      MFEM_ABORT("ADVectorFunction::operator(): This method should not be called. "
+                 "Use ADVectorFunction::operator(const ADVector &x, ADVector &F) instead.");
+   }
+   AD2Real_t operator()(const AD2Vector &x) const override final
+   {
+      MFEM_ABORT("ADVectorFunction::operator(): This method should not be called. "
+                 "Use ADVectorFunction::operator(const AD2Vector &x, AD2Vector &F) instead.");
+   }
+};
+
+class DifferentiableCoefficient : public Coefficient
+{
+private:
+   int idx; // index of the next input variable
+
+   class GradientCoefficient : public VectorCoefficient
+   {
+      DifferentiableCoefficient &c;
+   public:
+      GradientCoefficient(int dim, DifferentiableCoefficient &c)
+         : VectorCoefficient(dim), c(c) { }
+      void Eval(Vector &J, ElementTransformation &T,
+                const IntegrationPoint &ip) override
+      {
+         return c.f.Gradient(c.evaluator.Eval(T, ip), T, ip, J);
+      }
+   };
+
+   friend class GradientCoefficient;
+   GradientCoefficient grad_cf;
+
+   class HessianCoefficient : public MatrixCoefficient
+   {
+      DifferentiableCoefficient &c;
+   public:
+      HessianCoefficient(int dim, DifferentiableCoefficient &c)
+         : MatrixCoefficient(dim), c(c) { }
+      void Eval(DenseMatrix &H, ElementTransformation &T,
+                const IntegrationPoint &ip) override
+      { return c.f.Hessian(c.evaluator.Eval(T, ip), T, ip, H); }
+   };
+
+   friend class HessianCoefficient;
+   HessianCoefficient hess_cf;
+
+protected:
+   Evaluator evaluator;
+
+   ADFunction &f;
+public:
+   DifferentiableCoefficient(ADFunction &f)
+      : f(f), idx(0)
+      , grad_cf(f.n_input, *this)
+      , hess_cf(f.n_input, *this)
+   {}
+   DifferentiableCoefficient &AddInput(Evaluator::param_t param)
+   { evaluator.Add(param); return *this; }
+
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip) override
+   { return f(evaluator.Eval(T, ip), T, ip); }
+
+   GradientCoefficient& Gradient() { return grad_cf; }
+   HessianCoefficient& Hessian() { return hess_cf; }
+
+protected:
+};
+
+// Macro to generate type-varying implementation for ADFunction.
+// See, DiffusionEnergy, ..., for example of usage.
+// @param SCALAR is the name of templated scalar type
+// @param VEC is the name of templated vector type
+// @param MAT is the name of templated matrix type
+// @param var is the input variable name
+// @param body is the main function body. Use T() to create T-typed 0.
+#define AD_IMPL(SCALAR, VEC, MAT, var, body)                                           \
+   using ADFunction::operator();                                                       \
+   real_t operator()(const Vector &var) const override                                 \
+   {                                                                                   \
+      MFEM_ASSERT(var.Size() == n_input,                                               \
+                 "ADFunction::operator(): var.Size()=" << var.Size()                   \
+                  <<  " must match n_input=" << n_input)                               \
+      using SCALAR = real_t;                                                           \
+      using VEC = Vector;                                                              \
+      using MAT = DenseMatrix;                                                         \
+      body                                                                             \
+   }                                                                                   \
+                                                                                       \
+   ADReal_t operator()(const ADVector &var) const override                             \
+   {                                                                                   \
+      MFEM_ASSERT(var.Size() == n_input,                                               \
+                 "ADFunction::operator(): var.Size()=" << var.Size()                   \
+                  <<  " must match n_input=" << n_input)                               \
+      using SCALAR = ADReal_t;                                                         \
+      using VEC = ADVector;                                                            \
+      using MAT = ADMatrix;                                                            \
+      body                                                                             \
+   }                                                                                   \
+                                                                                       \
+   AD2Real_t operator()(const AD2Vector &var) const override                           \
+   {                                                                                   \
+      MFEM_ASSERT(var.Size() == n_input,                                               \
+                 "ADFunction::operator(): var.Size()=" << var.Size()                   \
+                  <<  " must match n_input=" << n_input)                               \
+      using SCALAR = AD2Real_t;                                                        \
+      using VEC = AD2Vector;                                                           \
+      using MAT = AD2Matrix;                                                           \
+      body                                                                             \
+   }
+
+
+// Macro to generate type-varying implementation for ADVectorFunction.
+// @param SCALAR is the name of templated scalar type
+// @param VEC is the name of templated vector type
+// @param MAT is the name of templated matrix type
+// @param var is the input variable name
+// @param result is the output variable name
+// @param body is the main function body. Use T() to create T-typed 0.
+#define AD_VEC_IMPL(SCALAR, VEC, MAT, var, result, body)                               \
+   using ADVectorFunction::operator();                                                 \
+   using ADVectorFunction::Gradient;                                                   \
+   using ADVectorFunction::Hessian;                                                    \
+                                                                                       \
+   void operator()(const Vector &var, Vector &result) const override                   \
+   {                                                                                   \
+      MFEM_ASSERT(var.Size() == n_input,                                               \
+                 "ADFunction::operator(): var.Size()=" << var.Size()                   \
+                  <<  " must match n_input=" << n_input)                               \
+      using SCALAR = real_t;                                                           \
+      using VEC = Vector;                                                              \
+      using MAT = DenseMatrix;                                                         \
+      body                                                                             \
+   }                                                                                   \
+                                                                                       \
+   void operator()(const ADVector &var, ADVector &result) const override               \
+   {                                                                                   \
+      MFEM_ASSERT(var.Size() == n_input,                                               \
+                 "ADFunction::operator(): var.Size()=" << var.Size()                   \
+                  <<  " must match n_input=" << n_input)                               \
+      using SCALAR = ADReal_t;                                                         \
+      using VEC = ADVector;                                                            \
+      using MAT = ADMatrix;                                                            \
+      body                                                                             \
+   }                                                                                   \
+                                                                                       \
+   void operator()(const AD2Vector &var, AD2Vector &result) const override             \
+   {                                                                                   \
+      MFEM_ASSERT(var.Size() == n_input,                                               \
+                 "ADFunction::operator(): var.Size()=" << var.Size()                   \
+                  <<  " must match n_input=" << n_input)                               \
+      using SCALAR = AD2Real_t;                                                        \
+      using VEC = AD2Vector;                                                           \
+      using MAT = AD2Matrix;                                                           \
+      body                                                                             \
+   }
+
+class MassEnergy : public ADFunction
+{
+public:
+   MassEnergy(int n_var)
+      : ADFunction(n_var)
+   {}
+   AD_IMPL(T, V, M, x, return 0.5*(x*x););
+};
+class DiffusionEnergy : public ADFunction
+{
+   const int dim;
+   mutable const Vector *K;
+public:
+   DiffusionEnergy(int dim)
+      : ADFunction(dim), dim(dim)
+   {}
+   DiffusionEnergy(int dim, Evaluator::param_t K)
+      : DiffusionEnergy(dim)
+   { SetK(K); }
+
+   void ProcessParameters(ElementTransformation &Tr,
+                          const IntegrationPoint &ip) const override
+   { K = &evaluator.Eval(Tr, ip); }
+
+   void SetK(Evaluator::param_t param)
+   {
+      int i = AddParameter(param);
+      int size = evaluator.val.GetBlock(i).Size();
+      MFEM_VERIFY(size == 1 || size == n_input || size == n_input*n_input,
+                  "Incorrect size for K. Dimension is " << n_input << "but K has size " << size);
+   }
+
+   AD_IMPL(T, V, M, gradu,
+   {
+      const int dim = gradu.Size();
+      const int Kdim = K->Size();
+      // No diffusion coefficient, ||grad u||^2
+      if (Kdim == 0) { return 0.5*(gradu*gradu); }
+      // Scalar diffusion coefficient, ||K^{1/2} grad u||^2
+      if (Kdim == 1) { return 0.5*(*K)[0]*(gradu*gradu); }
+      // Vector diffusion coefficient, ||diag(K)^{1/2} grad u||^2
+      if (Kdim == dim)
+      {
+         T result = T();
+         for (int i=0; i<dim; i++)
+         {
+            result += (*K)[i]*gradu[i]*gradu[i];
+         }
+         return 0.5*result;
+      }
+      // Matrix diffusion coefficient, ||K^{1/2} grad u||^2
+      if (Kdim == dim*dim)
+      {
+         DenseMatrix Kmat(K->GetData(), dim, dim);
+         T result = T();
+         for (int j=0; j<dim; j++)
+         {
+            for (int i=0; i<dim; i++)
+            {
+               result += Kmat(i,j)*gradu[i]*gradu[j];
+            }
+         }
+         return 0.5*result;
+      }
+      MFEM_ABORT("DiffusionEnergy: K must be a scalar, vector of size dim, "
+                 "or matrix of size dim x dim");
+      return T();
+   });
+};
+
+class DiffEnergy : public ADFunction
+{
+   const ADFunction &energy;
+   mutable const Vector *target;
+public:
+   DiffEnergy(const ADFunction &energy)
+      : ADFunction(energy.n_input)
+      , energy(energy)
+   { }
+
+   DiffEnergy(const ADFunction &energy, Evaluator::param_t other)
+      : DiffEnergy(energy)
+   {
+      int i = AddParameter(other);
+      MFEM_VERIFY(evaluator.val.GetBlock(0).Size() == n_input,
+                  "DiffEnergy: The provided target has the wrong size. "
+                  "Expected " << n_input << ", got " << evaluator.val.GetBlock(0).Size());
+   }
+
+   void SetTarget(Evaluator::param_t &target)
+   {
+      if (evaluator.val.NumBlocks() == 1)
+      { evaluator.Replace(0, target); }
+      else
+      { evaluator.Add(target); }
+      MFEM_VERIFY(evaluator.val.GetBlock(0).Size() == n_input,
+                  "DiffEnergy: The provided target has the wrong size. "
+                  "Expected " << n_input << ", got " << evaluator.val.GetBlock(0).Size());
+   }
+
+   void ProcessParameters(const BlockVector &x) const override
+   {
+      target = &x.GetBlock(0);
+   }
+
+   AD_IMPL(T, V, M, x,
+   {
+      V diff(x);
+      for (int i=0; i<n_input; i++)
+      { diff[i] -= (*target)[i]; }
+      return energy(diff);
+   });
+};
+
+class LinearElasticityEnergy : public ADFunction
+{
+   const int dim;
+   real_t &lambda;
+   real_t &mu;
+public:
+   void ProcessParameters(ElementTransformation &Tr,
+                          const IntegrationPoint &ip) const override
+   {
+      evaluator.Eval(Tr, ip);
+   }
+   LinearElasticityEnergy(int dim, Evaluator::param_t lambda,
+                          Evaluator::param_t mu, int offset=0)
+      : ADFunction(dim*dim, 2)
+      , dim(dim)
+      , lambda(*(evaluator.val.GetData() + offset))
+      , mu(*(evaluator.val.GetData() + evaluator.GetSize(mu) + offset))
+   {
+      int lambda_idx = evaluator.Add(lambda);
+      int mu_idx = evaluator.Add(mu);
+      MFEM_VERIFY(lambda_idx == 0,
+                  "LinearElasticityEnergy: lambda must be the first parameter");
+   }
+   AD_IMPL(T, V, M, gradu,
+   {
+      T divnorm = T();
+      for (int i=0; i<dim; i++) { divnorm += gradu[i*dim + i]; }
+      divnorm = divnorm*divnorm;
+      T h1_norm = T();
+      for (int i=0; i<dim; i++)
+      {
+         for (int j=0; j<dim; j++)
+         {
+            T symm = 0.5*(gradu[i*dim + j] + gradu[j*dim + i]);
+            h1_norm += symm*symm;
+         }
+      }
+      return 0.5*lambda*divnorm + mu*h1_norm;
+   });
+};
+
+// Lagrangian functional
+// f(x) + sum lambda[i]*c[i](x)
+class Lagrangian : public ADFunction
+{
+private:
+   enum { OBJONLY=-2, FULL=-1, CON=0};
+   int eval_mode =
+      FULL; // -2: objective, -1: full Lagrangian, >=0: constraint comp
+
+   ADFunction &objective; // f(x)
+
+   std::vector<ADFunction*> eq_con; // c[i](x)
+   Vector eq_rhs; // c[i](x) = con_target[i]
+public:
+
+   Lagrangian(ADFunction &objective, const int n_eq_con)
+      : ADFunction(objective.n_input+n_eq_con)
+      , objective(objective)
+   {}
+
+   Lagrangian AddEqConstraint(ADFunction &constraint,
+                              real_t target = 0.0);
+   Lagrangian SetEqRHS(int idx, real_t target) { eq_rhs[idx] = target; return *this; }
+
+   // return f(x) + sum lambda[i]*c[i](x)
+   void FullMode() { this->eval_mode = FULL; }
+   // return f(x)
+   void ObjectiveMode() { this->eval_mode = OBJONLY; }
+   // return c[i](x)
+   void EqConstraintMode(int comp)
+   {
+      MFEM_VERIFY(comp >= 0 && comp < eq_con.size(),
+                  "ALFunctional: comp must be in [0, n_input)");
+      this->eval_mode = comp;
+   }
+
+   void ProcessParameters(ElementTransformation &Tr,
+                          const IntegrationPoint &ip) const override;
+
+   AD_IMPL(T, V, M, x_and_lambda,
+   {
+      const V x(x_and_lambda.GetData(), objective.n_input);
+      const V lambda(x_and_lambda.GetData() + objective.n_input,
+                     eq_con.size());
+      if (eval_mode >= 0) { return (*eq_con[eval_mode])(x); }
+
+      T result = objective(x);
+      if (eval_mode == OBJONLY) { return result; } // only objective
+      for (int i=0; i<eq_con.size(); i++) { result += (*eq_con[i])(x)*lambda[i]; }
+      return result;
+   });
+
+private:
+};
+
+// Augmented Lagrangian functional
+class ALFunctional : public ADFunction
+{
+private:
+   enum { OBJONLY=-2, FULLAL=-1, CON=0};
+   int al_eval_mode = FULLAL; // -2: objective, -1: full AL, >=0: constraint comp
+
+   ADFunction &objective; // f(x)
+
+   std::vector<ADFunction*> eq_con; // c[i](x)
+   Vector eq_rhs; // c[i](x) = con_target[i]
+   Vector lambda; // Lagrange multipliers
+   real_t penalty=1.0; // penalty
+public:
+
+   ALFunctional(ADFunction &objective)
+      : ADFunction(objective.n_input)
+      , objective(objective)
+   {}
+
+   ALFunctional AddEqConstraint(ADFunction &constraint,
+                                real_t target = 0.0);
+   ALFunctional SetEqRHS(int idx, real_t target) { eq_rhs[idx] = target; return *this; }
+
+   void SetLambda(const Vector &lambda);
+   const Vector &GetLambda() const { return lambda; }
+   Vector &GetLambda() { return lambda; }
+
+   void SetPenalty(real_t mu);
+   real_t GetPenalty() const {return penalty; }
+   real_t &GetPenalty() { return penalty; }
+
+   // Full AL mode: f(x) + sum lambda[i]*c[i](x) + mu/2 * sum c[i](x)^2
+   void ALMode() { this->al_eval_mode = FULLAL; }
+   // Objective mode: f(x)
+   void ObjectiveMode() { this->al_eval_mode = OBJONLY; }
+   // Constraint mode: c[i](x)
+   void EqConstraintMode(int comp)
+   {
+      MFEM_VERIFY(comp >= 0 && comp < eq_con.size(),
+                  "ALFunctional: comp must be in [0, n_input)");
+      this->al_eval_mode = comp;
+   }
+
+   void ProcessParameters(ElementTransformation &Tr,
+                          const IntegrationPoint &ip) const override;
+
+   AD_IMPL(T, V, M, x,
+   {
+      if (al_eval_mode >= 0) { return evalAL<T>(x, al_eval_mode); }
+
+      T result = objective(x);
+      if (al_eval_mode == OBJONLY) { return result; } // only objective
+
+      for (int i=0; i<eq_con.size(); i++) { result += evalAL<T>(x, i); }
+
+      return result;
+   });
+
+private:
+   // Evaluate lambda*c(x) + (mu/2)*c(x)^2
+   template <typename T, typename V>
+   T evalAL(V &x, int idx) const
+   {
+      T cx = (*eq_con[idx])(x) - eq_rhs[idx];
+      if (al_eval_mode >= 0) { return cx; } // if non-negative, only c(x)
+      return cx*(lambda[idx] + penalty*0.5*cx);
+   }
+};
+// ------------------------------------------------------------------------------
+// Implement dual max/min
+// ------------------------------------------------------------------------------
+template <typename value_type, typename gradient_type, typename other_type>
+MFEM_HOST_DEVICE
+inline future::dual<value_type, gradient_type> max(
+   future::dual<value_type, gradient_type> a,
+   other_type b)
+{
+   if (a > b)
+   {
+      return a;
+   }
+   else if (a < b)
+   {
+      if constexpr (std::is_same<other_type, real_t>::value)
+      {
+         return future::dual<value_type, gradient_type> {b};
+      }
+      else
+      {
+         return b;
+      }
+   }
+   else
+   {
+      // If values are equal, return the average (subgradient)
+      return 0.5*(a + b);
+   }
+}
+
+template <typename value_type, typename gradient_type, typename other_type>
+MFEM_HOST_DEVICE
+inline future::dual<value_type, gradient_type> min(
+   future::dual<value_type, gradient_type> a,
+   other_type b)
+{
+   if (a < b)
+   {
+      return a;
+   }
+   else if (a > b)
+   {
+      if constexpr (std::is_same<other_type, real_t>::value)
+      {
+         return future::dual<value_type, gradient_type> {b};
+      }
+      else
+      {
+         return b;
+      }
+   }
+   else
+   {
+      // If values are equal, return the average (subgradient)
+      return 0.5*(a + b);
+   }
+}
+}

miniapps/autodiff/example.hpp

@@ -1,876 +0,0 @@
-// Copyright (c) 2010-2025, 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 ADEXAMPLE_HPP
-#define ADEXAMPLE_HPP
-
-#include "mfem.hpp"
-#include "admfem.hpp"
-#include <memory>
-#include <iostream>
-#include <fstream>
-
-namespace mfem
-{
-
-/// Example: Implementation of the residual evaluation for p-Laplacian
-/// problem. The residual is evaluated at the integration points for PDE
-/// parameters vparam and state fields (derivatives with respect to x,y,z and
-/// primal field) stored in vector uu.
-template<typename TDataType, typename TParamVector, typename TStateVector,
-         int residual_size, int state_size, int param_size>
-class MyResidualFunctor
-{
-public:
-   /// The operator returns the first derivative of the energy with respect to
-   /// all state variables. These are set in vector uu and consist of the
-   /// derivatives with respect to x,y,z and the primal field. The derivative is
-   /// stored in vector rr with length equal to the length of vector uu.
-   void operator()(TParamVector &vparam, TStateVector &uu, TStateVector &rr)
-   {
-      MFEM_ASSERT(residual_size==4,
-                  "PLaplacianResidual residual_size should be equal to 4!");
-      real_t pp = vparam[0];
-      real_t ee = vparam[1];
-      real_t ff = vparam[2];
-
-      // The vector rr holds the gradients of the following expression:
-      // (u_x^2+u_y^2+u_z^2+\varepsilon^2)^(p/2)-f.u,
-      // where u_x,u_y,u_z are the gradients of the scalar field u.
-      // The state vector is defined as uu=[u_x,u_y,u_z,u].
-
-      TDataType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
-      TDataType tvar = pow(ee * ee + norm2, (pp - 2.0) / 2.0);
-
-      rr[0] = tvar * uu[0];
-      rr[1] = tvar * uu[1];
-      rr[2] = tvar * uu[2];
-      rr[3] = -ff;
-   }
-};
-
-/// Defines template class (functor) for evaluating the energy of the
-/// p-Laplacian problem. The input parameters vparam are: vparam[0] - the
-/// p-Laplacian power, vparam[1] small value ensuring exciting of an unique
-/// solution, and vparam[2] - the distributed external input to the PDE. The
-/// template parameter TDataType will be replaced by the compiler with the
-/// appropriate AD type for automatic differentiation. The TParamVector
-/// represents the vector type used for the parameter vector, and TStateVector
-/// the vector type used for the state vector. The template parameters
-/// state_size and param_size provide information for the size of the state and
-/// the parameters vectors.
-template<typename TDataType, typename TParamVector, typename TStateVector
-         , int state_size, int param_size>
-class MyEnergyFunctor
-{
-public:
-   /// Returns the energy of a p-Laplacian for state field input provided in
-   /// vector uu and parameters provided in vector vparam.
-   TDataType operator()(TParamVector &vparam, TStateVector &uu)
-   {
-      MFEM_ASSERT(state_size==4,"MyEnergyFunctor state_size should be equal to 4!");
-      MFEM_ASSERT(param_size==3,"MyEnergyFunctor param_size should be equal to 3!");
-      real_t pp = vparam[0];
-      real_t ee = vparam[1];
-      real_t ff = vparam[2];
-
-      TDataType u = uu[3];
-      TDataType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
-
-      TDataType rez = pow(ee * ee + norm2, pp / 2.0) / pp - ff * u;
-      return rez;
-   }
-};
-
-
-/// Implements integrator for a p-Laplacian problem. The integrator is based on
-/// a class QFunction utilized for evaluating the energy, the first derivative
-/// (residual) and the Hessian of the energy (the Jacobian of the residual).
-/// The template parameter CQVectAutoDiff represents the automatically
-/// differentiated energy or residual implemented by the user.
-/// CQVectAutoDiff::VectorFunc(Vector parameters, Vector state,Vector residual)
-/// evaluates the residual at an integration point.
-/// CQVectAutoDiff::Jacobian(Vector parameters, Vector state, Matrix hessian)
-/// evaluates the Hessian of the energy(the Jacobian of the residual).
-template<class CQVectAutoDiff>
-class pLaplaceAD : public NonlinearFormIntegrator
-{
-protected:
-   Coefficient *pp;
-   Coefficient *coeff;
-   Coefficient *load;
-
-   CQVectAutoDiff rdf;
-
-public:
-   pLaplaceAD()
-   {
-      coeff = nullptr;
-      pp = nullptr;
-      load = nullptr;
-
-      vparam.SetSize(3);
-      vparam[0] = 2.0;  // default power
-      vparam[1] = 1e-8; // default epsilon
-      vparam[2] = 1.0;  // default load
-   }
-
-   pLaplaceAD(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr)
-   {
-      vparam.SetSize(3);
-      vparam[0] = 2.0;  // default power
-      vparam[1] = 1e-8; // default epsilon
-      vparam[2] = 1.0;  // default load
-
-   }
-
-   pLaplaceAD(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
-      : pp(&pp_), coeff(&q), load(&ld_)
-   {
-      vparam.SetSize(3);
-      vparam[0] = 2.0;  // default power
-      vparam[1] = 1e-8; // default epsilon
-      vparam[2] = 1.0;  // default load
-   }
-
-   virtual ~pLaplaceAD() {}
-
-   real_t GetElementEnergy(const FiniteElement &el,
-                           ElementTransformation &trans,
-                           const Vector &elfun) override
-   {
-      real_t energy = 0.0;
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      // derivatives in isoparametric coordinates
-      DenseMatrix dshape_iso(ndof, ndim);
-      // derivatives in physical space
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector grad(spaceDim);
-
-      Vector uu(4);     //[diff_x,diff_y,diff_z,u]
-
-      uu = 0.0;
-
-      // Calculates the functional/energy at an integration point.
-      MyEnergyFunctor<real_t,Vector,Vector,4,3> qfunc;
-
-      real_t w;
-      real_t detJ;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         detJ = (square ? w : w * w);
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         // AdjugateJacobian = / adj(J),         if J is square
-         //                    \ adj(J^t.J).J^t, otherwise
-         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-         // dshape_xyz should be divided by detJ for obtaining the real value
-         // calculate the gradient
-         dshape_xyz.MultTranspose(elfun, grad);
-
-         // set the power
-         if (pp != nullptr)
-         {
-            vparam[0] = pp->Eval(trans, ip);
-         }
-
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            vparam[1] = coeff->Eval(trans, ip);
-         }
-         // add the contribution from the load
-         if (load != nullptr)
-         {
-            vparam[2] = load->Eval(trans, ip);
-         }
-         // fill the values of vector uu
-         for (int jj = 0; jj < spaceDim; jj++)
-         {
-            uu[jj] = grad[jj] / detJ;
-         }
-         uu[3] = shapef * elfun;
-         // the energy is taken directly from the templated function
-         energy = energy + w * qfunc(vparam,uu);
-      }
-      return energy;
-   }
-
-   void AssembleElementVector(const FiniteElement &el,
-                              ElementTransformation &trans,
-                              const Vector &elfun,
-                              Vector &elvect) override
-   {
-      MFEM_PERF_BEGIN("AssembleElementVector");
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector lvec(ndof);
-      elvect.SetSize(ndof);
-      elvect = 0.0;
-
-      DenseMatrix B(ndof, 4); // [diff_x,diff_y,diff_z, shape]
-      Vector uu(4);           // [diff_x,diff_y,diff_z,u]
-      Vector du(4);
-      B = 0.0;
-      uu = 0.0;
-      real_t w;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
-
-         // set the matrix B
-         for (int jj = 0; jj < spaceDim; jj++)
-         {
-            B.SetCol(jj, dshape_xyz.GetColumn(jj));
-         }
-         B.SetCol(3, shapef);
-
-         // set the power
-         if (pp != nullptr)
-         {
-            vparam[0] = pp->Eval(trans, ip);
-         }
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            vparam[1] = coeff->Eval(trans, ip);
-         }
-         // add the contribution from the load
-         if (load != nullptr)
-         {
-            vparam[2] = load->Eval(trans, ip);
-         }
-
-         // calculate uu
-         B.MultTranspose(elfun, uu);
-         // calculate derivative of the energy with respect to uu
-         rdf.VectorFunc(vparam,uu,du);
-         B.Mult(du, lvec);
-         elvect.Add(w, lvec);
-      } // end integration loop
-      MFEM_PERF_END("AssembleElementVector");
-   }
-
-   void AssembleElementGrad(const FiniteElement &el,
-                            ElementTransformation &trans,
-                            const Vector &elfun,
-                            DenseMatrix &elmat) override
-   {
-      MFEM_PERF_BEGIN("AssembleElementGrad");
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      elmat.SetSize(ndof, ndof);
-      elmat = 0.0;
-
-      DenseMatrix B(ndof, 4); // [diff_x,diff_y,diff_z, shape]
-      DenseMatrix A(ndof, 4);
-      Vector uu(4); // [diff_x,diff_y,diff_z,u]
-      DenseMatrix duu(4, 4);
-      B = 0.0;
-      uu = 0.0;
-
-      real_t w;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
-
-         // set the matrix B
-         for (int jj = 0; jj < spaceDim; jj++)
-         {
-            B.SetCol(jj, dshape_xyz.GetColumn(jj));
-         }
-         B.SetCol(3, shapef);
-
-         // set the power
-         if (pp != nullptr)
-         {
-            vparam[0] = pp->Eval(trans, ip);
-         }
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            vparam[1] = coeff->Eval(trans, ip);
-         }
-         // add the contribution from the load
-         if (load != nullptr)
-         {
-            vparam[2] = load->Eval(trans, ip);
-         }
-
-         // calculate uu
-         B.MultTranspose(elfun, uu);
-         // calculate derivative of the energy with respect to uu
-         rdf.Jacobian(vparam,uu,duu);
-         Mult(B, duu, A);
-         AddMult_a_ABt(w, A, B, elmat);
-
-      } // end integration loop
-      MFEM_PERF_END("AssembleElementGrad");
-   }
-
-private:
-   Vector vparam; // [power, epsilon, load]
-
-};
-
-/// Implements hand-coded integrator for a p-Laplacian problem. Utilized as
-/// alternative for the pLaplaceAD class based on automatic differentiation.
-class pLaplace : public NonlinearFormIntegrator
-{
-protected:
-   Coefficient *pp;
-   Coefficient *coeff;
-   Coefficient *load;
-
-public:
-   pLaplace()
-   {
-      coeff = nullptr;
-      pp = nullptr;
-      load = nullptr;
-   }
-
-   pLaplace(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
-
-   pLaplace(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
-      : pp(&pp_), coeff(&q), load(&ld_)
-   {}
-
-   virtual ~pLaplace() {}
-
-   real_t GetElementEnergy(const FiniteElement &el,
-                           ElementTransformation &trans,
-                           const Vector &elfun) override
-   {
-      real_t energy = 0.0;
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector grad(spaceDim);
-
-      real_t w;
-      real_t detJ;
-      real_t nrgrad2;
-      real_t ppp = 2.0;
-      real_t eee = 0.0;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         detJ = (square ? w : w * w);
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         // AdjugateJacobian = / adj(J),         if J is square
-         //                    \ adj(J^t.J).J^t, otherwise
-         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-         // dshape_xyz should be divided by detJ for obtaining the real value
-         // calculate the gradient
-         dshape_xyz.MultTranspose(elfun, grad);
-         nrgrad2 = grad * grad / (detJ * detJ);
-
-         // set the power
-         if (pp != nullptr)
-         {
-            ppp = pp->Eval(trans, ip);
-         }
-
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            eee = coeff->Eval(trans, ip);
-         }
-
-         energy = energy + w * std::pow(nrgrad2 + eee * eee, ppp / 2.0) / ppp;
-
-         // add the contribution from the load
-         if (load != nullptr)
-         {
-            energy = energy - w * (shapef * elfun) * load->Eval(trans, ip);
-         }
-      }
-      return energy;
-   }
-
-   void AssembleElementVector(const FiniteElement &el,
-                              ElementTransformation &trans,
-                              const Vector &elfun,
-                              Vector &elvect) override
-   {
-      MFEM_PERF_BEGIN("AssembleElementVector");
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector grad(spaceDim);
-      Vector lvec(ndof);
-      elvect.SetSize(ndof);
-      elvect = 0.0;
-
-      real_t w;
-      real_t detJ;
-      real_t nrgrad;
-      real_t aa;
-      real_t ppp = 2.0;
-      real_t eee = 0.0;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         detJ = (square ? w : w * w);
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         // AdjugateJacobian = / adj(J),         if J is square
-         //                    \ adj(J^t.J).J^t, otherwise
-         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-         // dshape_xyz should be divided by detJ for obtaining the real value
-
-         // calculate the gradient
-         dshape_xyz.MultTranspose(elfun, grad);
-         nrgrad = grad.Norml2() / detJ;
-         // grad is not scaled so far, i.e., grad=grad/detJ
-
-         // set the power
-         if (pp != nullptr)
-         {
-            ppp = pp->Eval(trans, ip);
-         }
-
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            eee = coeff->Eval(trans, ip);
-         }
-         // compute (norm of the gradient)^2 + epsilon^2
-         aa = nrgrad * nrgrad + eee * eee;
-         aa = std::pow(aa, (ppp - 2.0) / 2.0);
-         dshape_xyz.Mult(grad, lvec);
-         elvect.Add(w * aa / (detJ * detJ), lvec);
-
-         // add loading
-         if (load != nullptr)
-         {
-            elvect.Add(-w * load->Eval(trans, ip), shapef);
-         }
-      } // end integration loop
-      MFEM_PERF_END("AssembleElementVector");
-   }
-
-   void AssembleElementGrad(const FiniteElement &el,
-                            ElementTransformation &trans,
-                            const Vector &elfun,
-                            DenseMatrix &elmat) override
-   {
-      MFEM_PERF_BEGIN("AssembleElementGrad");
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector grad(spaceDim);
-      Vector lvec(ndof);
-      // set the size of the element matrix
-      elmat.SetSize(ndof, ndof);
-      elmat = 0.0;
-
-      real_t w; // integration weight
-      real_t detJ;
-      real_t nrgrad; // norm of the gradient
-      real_t aa0; // original nonlinear diffusion coefficient
-      real_t aa1; // gradient of the above
-      real_t ppp = 2.0; // power in the P-Laplacian
-      real_t eee = 0.0; // regularization parameter
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         detJ = (square ? w : w * w);
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         // AdjugateJacobian = / adj(J),         if J is square
-         //                    \ adj(J^t.J).J^t, otherwise
-         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-         // dshape_xyz should be divided by detJ for obtaining the real value
-         // grad is not scaled so far,i.e., grad=grad/detJ
-
-         // set the power
-         if (pp != nullptr)
-         {
-            ppp = pp->Eval(trans, ip);
-         }
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            eee = coeff->Eval(trans, ip);
-         }
-
-         // calculate the gradient
-         dshape_xyz.MultTranspose(elfun, grad);
-         nrgrad = grad.Norml2() / detJ;
-         // (u_x^2+u_y^2+u_z^2+\varepsilon^2)
-         aa0 = nrgrad * nrgrad + eee * eee;
-         aa1 = std::pow(aa0, (ppp - 2.0) / 2.0);
-         aa0 = (ppp - 2.0) * std::pow(aa0, (ppp - 4.0) / 2.0);
-         dshape_xyz.Mult(grad, lvec);
-         w = w / (detJ * detJ);
-         AddMult_a_VVt(w * aa0 / (detJ * detJ), lvec, elmat);
-         AddMult_a_AAt(w * aa1, dshape_xyz, elmat);
-
-      } // end integration loop
-      MFEM_PERF_END("AssembleElementGrad");
-   }
-};
-
-/// Implements AD enabled integrator for a p-Laplacian problem. The tangent
-/// matrix is computed using the residual of the element. The template argument
-/// should be equal to the size of the residual vector (element vector), i.e.,
-/// the user should specify the size to match the exact vector size for the
-/// considered order of the shape functions.
-
-template<int sizeres=10>
-class pLaplaceSL : public NonlinearFormIntegrator
-{
-protected:
-   Coefficient *pp;
-   Coefficient *coeff;
-   Coefficient *load;
-
-public:
-   pLaplaceSL()
-   {
-      coeff = nullptr;
-      pp = nullptr;
-      load = nullptr;
-   }
-
-   pLaplaceSL(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
-
-   pLaplaceSL(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
-      : pp(&pp_), coeff(&q), load(&ld_)
-   {}
-
-   virtual ~pLaplaceSL() {}
-
-   real_t GetElementEnergy(const FiniteElement &el,
-                           ElementTransformation &trans,
-                           const Vector &elfun) override
-   {
-      real_t energy = 0.0;
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector grad(spaceDim);
-
-      real_t w;
-      real_t detJ;
-      real_t nrgrad2;
-      real_t ppp = 2.0;
-      real_t eee = 0.0;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         detJ = (square ? w : w * w);
-         w = ip.weight * w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         // AdjugateJacobian = / adj(J),         if J is square
-         //                    \ adj(J^t.J).J^t, otherwise
-         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-         // dshape_xyz should be divided by detJ for obtaining the real value
-         // calculate the gradient
-         dshape_xyz.MultTranspose(elfun, grad);
-         nrgrad2 = grad * grad / (detJ * detJ);
-
-         // set the power
-         if (pp != nullptr)
-         {
-            ppp = pp->Eval(trans, ip);
-         }
-
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            eee = coeff->Eval(trans, ip);
-         }
-
-         energy = energy + w * std::pow(nrgrad2 + eee * eee, ppp / 2.0) / ppp;
-
-         // add the contribution from the load
-         if (load != nullptr)
-         {
-            energy = energy - w * (shapef * elfun) * load->Eval(trans, ip);
-         }
-      }
-      return energy;
-   }
-
-   void AssembleElementVector(const FiniteElement &el,
-                              ElementTransformation &trans,
-                              const Vector &elfun,
-                              Vector &elvect) override
-   {
-      MFEM_PERF_BEGIN("AssembleElementVector");
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      Vector shapef(ndof);
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      Vector grad(spaceDim);
-      Vector lvec(ndof);
-      elvect.SetSize(ndof);
-      elvect = 0.0;
-
-      real_t w;
-      real_t detJ;
-      real_t nrgrad;
-      real_t aa;
-      real_t ppp = 2.0;
-      real_t eee = 0.0;
-
-      for (int i = 0; i < ir.GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(i);
-         trans.SetIntPoint(&ip);
-         w = trans.Weight();
-         detJ = (square ? w : w * w);
-         w = ip.weight * w; //w;
-
-         el.CalcDShape(ip, dshape_iso);
-         el.CalcShape(ip, shapef);
-         // AdjugateJacobian = / adj(J),         if J is square
-         //                    \ adj(J^t.J).J^t, otherwise
-         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-         // dshape_xyz should be divided by detJ for obtaining the real value
-
-         // calculate the gradient
-         dshape_xyz.MultTranspose(elfun, grad);
-         nrgrad = grad.Norml2() / detJ;
-         // grad is not scaled so far, i.e., grad=grad/detJ
-
-         // set the power
-         if (pp != nullptr)
-         {
-            ppp = pp->Eval(trans, ip);
-         }
-
-         // set the coefficient ensuring positiveness of the tangent matrix
-         if (coeff != nullptr)
-         {
-            eee = coeff->Eval(trans, ip);
-         }
-
-         aa = nrgrad * nrgrad + eee * eee;
-         aa = std::pow(aa, (ppp - 2.0) / 2.0);
-         dshape_xyz.Mult(grad, lvec);
-         elvect.Add(w * aa / (detJ * detJ), lvec);
-
-         // add loading
-         if (load != nullptr)
-         {
-            elvect.Add(-w * load->Eval(trans, ip), shapef);
-         }
-      } // end integration loop
-      MFEM_PERF_END("AssembleElementVector");
-   }
-
-   void AssembleElementGrad(const FiniteElement &el,
-                            ElementTransformation &trans,
-                            const Vector &elfun,
-                            DenseMatrix &elmat) override
-   {
-      MFEM_PERF_BEGIN("AssembleElementGrad");
-      const int ndof = el.GetDof();
-      const int ndim = el.GetDim();
-      const int spaceDim = trans.GetSpaceDim();
-      bool square = (ndim == spaceDim);
-      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
-      const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
-
-      DenseMatrix dshape_iso(ndof, ndim);
-      DenseMatrix dshape_xyz(ndof, spaceDim);
-      elmat.SetSize(ndof, ndof);
-      elmat = 0.0;
-
-      real_t w;
-      real_t detJ;
-      real_t ppp = 2.0;
-      real_t eee = 0.0;
-
-      mfem::Vector param(3); param=0.0;
-
-      // Computes the residual at an integration point. The implementation is a
-      // copy of the integration loop in AssembleElementVector.
-      auto resfun = [&](mfem::Vector& vparam, mfem::ad::ADVectorType& uu,
-                        mfem::ad::ADVectorType& vres)
-      {
-
-         vres.SetSize(uu.Size()); vres=0.0;
-         mfem::ad::ADVectorType grad(spaceDim);
-         mfem::ad::ADFloatType nrgrad;
-         mfem::ad::ADFloatType aa;
-         mfem::ad::ADVectorType lvec(ndof);
-
-         for (int q = 0; q < ir.GetNPoints(); q++)
-         {
-            lvec=0.0;
-
-            const IntegrationPoint &ip = ir.IntPoint(q);
-            trans.SetIntPoint(&ip);
-            w = trans.Weight();
-            detJ = (square ? w : w * w);
-            w = ip.weight * w;
-
-            el.CalcDShape(ip, dshape_iso);
-            // AdjugateJacobian = / adj(J),         if J is square
-            //                    \ adj(J^t.J).J^t, otherwise
-            Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
-            // dshape_xyz should be divided by detJ for obtaining the real value
-            // grad is not scaled so far,i.e., grad=grad/detJ
-
-            // set the power
-            if (pp != nullptr)
-            {
-               ppp = pp->Eval(trans, ip);
-            }
-            // set the coefficient ensuring positiveness of the tangent matrix
-            if (coeff != nullptr)
-            {
-               eee = coeff->Eval(trans, ip);
-            }
-
-            grad=0.0;
-            // calculate the gradient
-            for (int i=0; i<spaceDim; i++)
-            {
-               for (int j=0; j<ndof; j++)
-               {
-                  grad[i]= grad[i]+ dshape_xyz(j,i)*uu[j];
-               }
-            }
-
-
-            nrgrad= (grad*grad)/(detJ*detJ);
-
-            aa = nrgrad + eee * eee;
-            aa = pow(aa, (ppp - 2.0) / 2.0);
-
-            for (int i=0; i<spaceDim; i++)
-            {
-               for (int j=0; j<ndof; j++)
-               {
-                  lvec[j] = lvec[j] + dshape_xyz(j,i) * grad[i];
-               }
-            }
-
-            for (int j=0; j<ndof; j++)
-            {
-               vres[j]=vres[j] + lvec[j] * (w*aa/(detJ*detJ));
-            }
-         }
-      };
-
-      mfem::Vector bla(elfun);
-      // calculate the gradient - only for a fixed ndof
-      mfem::VectorFuncAutoDiff<sizeres,sizeres,3> fdr(resfun);
-      fdr.Jacobian(param, bla, elmat);
-      MFEM_PERF_END("AssembleElementGrad");
-   }
-};
-
-} // namespace mfem
-
-#endif

miniapps/autodiff/logger.cpp

@@ -0,0 +1,261 @@
+#include "logger.hpp"
+
+namespace mfem
+{
+
+TableLogger::TableLogger(std::ostream &os)
+   : os(os), w(14), var_name_printed(false)
+{
+#ifdef MFEM_USE_MPI
+   isRoot = mfem::Mpi::IsInitialized() ? mfem::Mpi::Root() : true;
+#endif
+}
+
+void TableLogger::Append(const std::string name, double &val)
+{
+   names.push_back(name);
+   data_double.push_back(&val);
+   data_order.push_back(dtype::DOUBLE);
+}
+
+void TableLogger::Append(const std::string name, int &val)
+{
+   names.push_back(name);
+   data_int.push_back(&val);
+   data_order.push_back(dtype::INT);
+}
+
+void TableLogger::Print(bool print_varname)
+{
+   if (isRoot)
+   {
+      if (!var_name_printed || print_varname)
+      {
+         for (auto &name : names)
+         {
+            os << std::setw(w) << std::setfill(' ') << name << ",\t";
+         }
+         os << "\b\b";
+         os << std::endl;
+         if (!var_name_printed && file && file->is_open())
+         {
+            for (int i=0; i<names.size() - 1; i++)
+            {
+               *file << std::setw(w) << std::setfill(' ') << names[i] << ",\t";
+            }
+            *file << std::setw(w) << std::setfill(' ') << names.back() << std::endl;
+         }
+         var_name_printed = true;
+      }
+      int i(0), i_double(0), i_int(0);
+      for (int i=0; i<data_order.size(); i++)
+      {
+         auto d = data_order[i];
+         switch (d)
+         {
+            case dtype::DOUBLE:
+            {
+               os << std::setw(w) << *data_double[i_double];
+               if (file && file->is_open())
+               {
+                  *file << std::setprecision(8) << std::scientific << std::setw(w)
+                        << std::setfill(' ') << *data_double[i_double];
+               }
+               i_double++;
+               break;
+            }
+            case dtype::INT:
+            {
+               os << std::setw(w) << *data_int[i_int];
+               if (file && file->is_open())
+               {
+                  *file << std::setw(w) << std::setfill(' ') << *data_int[i_int];
+               }
+               i_int++;
+               break;
+            }
+            default:
+            {
+               MFEM_ABORT("Unknown data type. See, TableLogger::dtype");
+            }
+         }
+         if (i < data_order.size() - 1)
+         {
+            os << ",\t";
+            *file << ",\t";
+         }
+      }
+      os << std::endl;
+      if (file)
+      {
+         *file << std::endl;
+      }
+   }
+}
+
+void TableLogger::SaveWhenPrint(std::string filename, std::ios::openmode mode)
+{
+   if (isRoot)
+   {
+      filename = filename.append(".csv");
+      file.reset(new std::fstream);
+      file->open(filename, mode);
+      if (!file->is_open())
+      {
+         std::string msg("");
+         msg += "Cannot open file ";
+         msg += filename;
+         MFEM_ABORT(msg);
+      }
+   }
+}
+
+bool GLVis::Append(GridFunction *gf, QuadratureFunction *qf,
+                   std::string_view window_title, std::string_view keys)
+{
+   MFEM_VERIFY((gf == nullptr && qf != nullptr)
+               || (gf != nullptr && qf == nullptr),
+               "Either GridFunction or QuadratureFunction must be provided, "
+               "but not both.");
+   bool is_gf = gf != nullptr;
+   sockets.push_back(std::make_unique<socketstream>(hostname, port, secure));
+   socketstream &socket = *sockets.back();
+   if (!socket.is_open() || !socket.good())
+   {
+      MFEM_WARNING("GLVis: Cannot connect to " << hostname << ":" << port);
+      sockets.back().reset();
+      sockets.pop_back();
+      return false;
+   }
+   socket.precision(8);
+   gfs.Append(gf);
+   qfs.Append(qf);
+
+   Mesh *mesh;
+   if (is_gf) { mesh = gf->FESpace()->GetMesh(); }
+   else { mesh = qf->GetSpace()->GetMesh(); }
+   meshes.Append(mesh);
+
+   cfs.Append(nullptr);
+   vcfs.Append(nullptr);
+
+#ifdef MFEM_USE_MPI
+   parallel.Append(false);
+   myrank.Append(0);
+   nrrank.Append(1);
+   ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
+   if (pmesh != nullptr)
+   {
+      parallel.Last() = true;
+      nrrank.Last() = pmesh->GetNRanks();
+      myrank.Last() = pmesh->GetMyRank();
+      socket << "parallel " << nrrank.Last() << " " << myrank.Last() <<
+             "\n";
+   }
+#endif
+   if (is_gf)
+   {
+      socket << "solution\n" << *mesh << *gf;
+   }
+   else
+   {
+      socket << "quadrature\n" << *mesh << *qf << "\n";
+   }
+
+
+   if (!keys.empty())
+   {
+      socket << "keys " << keys << "\n";
+      bool hasQ=false;
+      if (!is_gf)
+      {
+         auto end_pos = std::min(keys.find(' '), keys.find('\n'));
+         std::string_view actual_keys = keys.substr(0, end_pos);
+         if (actual_keys.find('Q') != std::string_view::npos) { hasQ = true; }
+      }
+      qfkey_has_Q.Append(hasQ);
+   }
+   if (!window_title.empty())
+   {
+      socket << "window_title '" << window_title <<"'\n";
+   }
+   int row = (sockets.size() - 1) / nrWinPerRow;
+   int col = (sockets.size() - 1) % nrWinPerRow;
+   socket << " window_geometry "
+          << w*col << " " << h*row << " "
+          << w << " " << h << "\n";
+   socket << std::flush;
+#ifdef MFEM_USE_MPI
+   if (parallel.Last())
+   {
+      MPI_Comm comm = static_cast<ParMesh*>(meshes.Last())->GetComm();
+      MPI_Barrier(comm);
+   }
+#endif
+   return true;
+}
+
+void GLVis::Append(Coefficient &cf, QuadratureSpace &qs,
+                   std::string_view window_title,
+                   std::string_view keys)
+{
+   owned_qfs.push_back(std::make_unique<QuadratureFunction>(qs));
+   cf.Project(*owned_qfs.back());
+   if (Append(nullptr, owned_qfs.back().get(), window_title, keys))
+   {
+      cfs.Last() = &cf;
+   }
+}
+
+void GLVis::Append(VectorCoefficient &cf, QuadratureSpace &qs,
+                   std::string_view window_title,
+                   std::string_view keys)
+{
+   owned_qfs.push_back(std::make_unique<QuadratureFunction>(qs, cf.GetVDim()));
+   cf.Project(*owned_qfs.back());
+   if (Append(nullptr, owned_qfs.back().get(), window_title, keys))
+   {
+      vcfs.Last() = &cf;
+   }
+}
+
+void GLVis::Update()
+{
+   for (int i=0; i<sockets.size(); i++)
+   {
+      if (!sockets[i]->is_open() || !sockets[i]->good())
+      {
+         MFEM_WARNING("GLVis: Connection to " << hostname << ":" << port
+                      << " for window " << i+1 << " lost.");
+         continue;
+      }
+#ifdef MFEM_USE_MPI
+      if (parallel[i])
+      {
+         *sockets[i] << "parallel " << nrrank[i] << " " << myrank[i] <<
+                     "\n";
+      }
+#endif
+      if (gfs[i] != nullptr)
+      {
+         *sockets[i] << "solution\n" << *meshes[i] << *gfs[i];
+      }
+      else if (qfs[i] != nullptr)
+      {
+         if (cfs[i] != nullptr) { cfs[i]->Project(*qfs[i]); }
+         else if (vcfs[i] != nullptr) { vcfs[i]->Project(*qfs[i]); }
+         *sockets[i] << "quadrature\n" << *meshes[i] << *qfs[i];
+         if (qfkey_has_Q[i]) { *sockets[i] << "keys QQQ\n"; }
+      }
+      *sockets[i] << std::flush;
+#ifdef MFEM_USE_MPI
+      if (parallel[i])
+      {
+         MPI_Comm comm = static_cast<ParMesh*>(meshes[i])->GetComm();
+         MPI_Barrier(comm);
+      }
+#endif
+   }
+}
+
+} // namespace mfem

miniapps/autodiff/logger.hpp

@@ -0,0 +1,122 @@
+#pragma once
+
+#include <iomanip>
+#include <iostream>
+#include <vector>
+#include "mfem.hpp"
+
+
+namespace mfem
+{
+
+class TableLogger
+{
+public:
+   enum dtype { DOUBLE, INT };
+
+protected:
+   // Double data to be printed.
+   std::vector<double *> data_double;
+   // Int data to be printed
+   std::vector<int *> data_int;
+   // Data type for each column
+   std::vector<dtype> data_order;
+   // Name of each monitored data
+   std::vector<std::string> names;
+   // Output stream
+   std::ostream &os;
+   // Column width
+   int w;
+   // Whether the variable name row has been printed or not
+   bool var_name_printed;
+#ifdef MFEM_USE_MPI
+   bool isRoot; // true if serial or root in parallel
+#else
+   static constexpr bool isRoot = true;
+#endif
+   std::unique_ptr<std::fstream> file;
+
+private:
+public:
+   // Create a logger that prints a row of variables for each call of Print
+   TableLogger(std::ostream &os = std::cout);
+   // Set column width of the table to be printed
+   void setw(const int column_width) { w = column_width; }
+   // Add double data to be monitored
+   void Append(const std::string name, double &val);
+   // Add double data to be monitored
+   void Append(const std::string name, int &val);
+   // Print a row of currently monitored data. If it is called
+   void Print(bool print_valname=false);
+   // Save data to a file whenever Print is called.
+   void SaveWhenPrint(std::string filename,
+                      std::ios::openmode mode = std::ios::out);
+   // Close file manually.
+   void CloseFile() { if (file) { file.reset(nullptr); } }
+};
+
+class GLVis
+{
+   std::vector<std::unique_ptr<socketstream>> sockets;
+   // Array<mfem::socketstream *> sockets;
+   Array<mfem::GridFunction *> gfs;
+   Array<mfem::QuadratureFunction *> qfs;
+   Array<bool> qfkey_has_Q;
+   Array<bool> qfhas_cf;
+   Array<Coefficient*> cfs;
+   Array<VectorCoefficient*> vcfs;
+   std::vector<std::unique_ptr<QuadratureFunction>> owned_qfs;
+   Array<Mesh *> meshes;
+   Array<bool> parallel;
+   Array<int> myrank;
+   Array<int> nrrank;
+   const char *hostname;
+   const int port;
+   int w, h, nrWinPerRow;
+   bool secure;
+   bool Append(GridFunction *gf, QuadratureFunction *qf,
+               std::string_view window_title, std::string_view keys);
+
+public:
+#ifdef MFEM_USE_GNUTLS
+   static const bool secure_default = true;
+#else
+   static const bool secure_default = false;
+#endif
+   GLVis(const char hostname[], int port, int w=400, int h=350,
+         int nrWinPerRow=1,
+         bool secure = secure_default)
+      : sockets(0), gfs(0), meshes(0), parallel(0), hostname(hostname),
+        port(port), w(w), h(h), nrWinPerRow(nrWinPerRow),
+        secure(secure_default) {}
+
+   void Append(GridFunction &gf,
+               std::string_view window_title= {},
+               std::string_view keys= {})
+   { Append(&gf, nullptr, window_title, keys); }
+   void Append(QuadratureFunction &qf,
+               std::string_view window_title= {},
+               std::string_view keys= {})
+   { Append(nullptr, &qf, window_title, keys); }
+   void Append(Coefficient &cf, QuadratureSpace &qs,
+               std::string_view window_title= {},
+               std::string_view keys= {});
+   void Append(VectorCoefficient &cf, QuadratureSpace &qs,
+               std::string_view window_title= {},
+               std::string_view keys= {});
+   void Update();
+
+   GridFunction& GetGridFunction(int i)
+   {
+      MFEM_VERIFY(i < gfs.Size(), "Index out of range");
+      return *gfs[i];
+   }
+
+   socketstream &GetSocket(int i)
+   {
+      MFEM_VERIFY(i < sockets.size(), "Index out of range");
+      return *sockets[i];
+   }
+};
+
+} // namespace mfem

miniapps/autodiff/makefile

@@ -13,7 +13,7 @@
 MFEM_DIR ?= ../..
 MFEM_BUILD_DIR ?= ../..
 MFEM_INSTALL_DIR ?= ../../mfem
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/autodiff/,)
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/electromagnetics/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
@@ -21,49 +21,41 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
 include $(DEFAULTS_MK)
 
+AD_COMMON_SRC = logger.cpp ad_native.cpp pg.cpp
+AD_COMMON_OBJ = $(AD_COMMON_SRC:.cpp=.o)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-ADIFF_COMMON_SRC =
-ADIFF_COMMON_OBJ = $(ADIFF_COMMON_SRC:.cpp=.o)
+SEQ_MINIAPPS = ad_ex0 ad_ex1 ad_ex2 ad_ex3
+PAR_MINIAPPS = ad_ex6
+ifeq ($(MFEM_USE_PETSC),YES)
+		PAR_MINIAPPS += ad_ex4 ad_ex5
+endif
 
-SEQ_MINIAPPS = seq_example seq_test
-PAR_MINIAPPS = par_example
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
-   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+   MINIAPPS = $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 endif
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
 
 # Remove built-in rules
 %: %.cpp
 %.o: %.cpp
 
-%: %.o $(ADIFF_COMMON_OBJ)
+%: %.o $(AD_COMMON_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $^ -o $@ $(MFEM_LIBS)
-
 %.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+	$(MFEM_CXX) $(MFEM_FLAGS) -I$(MFEM_DIR)/miniapps/autodiff -c $< -o $@
 
 all: $(MINIAPPS)
 
-MFEM_TESTS = MINIAPPS
-include $(MFEM_TEST_MK)
-
-# Testing: Parallel vs. serial runs
-RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-TEST_NAME := ADIFF miniapp
-%-test-par: %
-	@$(call mfem-test,$<, $(RUN_MPI), $(TEST_NAME))
-%-test-seq: %
-	@$(call mfem-test,$<,, $(TEST_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)
@@ -72,7 +64,5 @@ clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf ParaView
 	rm -rf *.dSYM *.TVD.*breakpoints
-
-clean-exec:
-	@rm -rf Example* 

miniapps/autodiff/par_example.cpp

@@ -1,550 +0,0 @@
-// Copyright (c) 2010-2025, 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.
-//
-//                       MFEM AD Example  - Parallel Version
-//
-// Compile with: make par_example
-//
-// Sample runs:  mpirun -np 2 par_example -m ../../data/beam-quad.mesh -pp 3.8
-//               mpirun -np 2 par_example -m ../../data/beam-tri.mesh  -pp 7.2
-//               mpirun -np 2 par_example -m ../../data/beam-hex.mesh
-//               mpirun -np 2 par_example -m ../../data/beam-tet.mesh
-//               mpirun -np 2 par_example -m ../../data/beam-wedge.mesh
-//
-// Description:  This examples solves a quasi-static nonlinear p-Laplacian
-//               problem with zero Dirichlet boundary conditions applied on all
-//               defined boundaries
-//
-//               The example demonstrates the use of nonlinear operators
-//               combined with automatic differentiation (AD). The integrators
-//               are defined in example.hpp. Selecting integrator = 0 will use
-//               the manually implemented integrator. Selecting integrator = 1
-//               or 2 will utilize one of the AD integrators.
-//
-//               We recommend viewing examples 1 and 19, before viewing this
-//               example.
-
-#include "example.hpp"
-
-using namespace mfem;
-
-enum IntegratorType
-{
-   HandCodedIntegrator  = 0,
-   ADJacobianIntegrator = 1,
-   ADHessianIntegrator  = 2
-};
-
-/// Non-linear solver for the p-Laplacian problem.
-class ParNLSolverPLaplacian
-{
-public:
-   /// Constructor Input: imesh - FE mesh, finite element space, power for the
-   /// p-Laplacian, external load (source, input), regularization parameter
-   ParNLSolverPLaplacian(MPI_Comm comm, ParMesh& imesh,
-                         ParFiniteElementSpace& ifespace,
-                         real_t powerp=2,
-                         Coefficient* load=nullptr,
-                         real_t regularizationp=1e-7)
-   {
-      lcomm = comm;
-
-      // default parameters for the Newton solver
-      newton_rtol = 1e-4;
-      newton_atol = 1e-8;
-      newton_iter = 10;
-
-      // linear solver
-      linear_rtol = 1e-7;
-      linear_atol = 1e-15;
-      linear_iter = 500;
-
-      print_level = 0;
-
-      // set the mesh
-      mesh=&imesh;
-
-      // set the fespace
-      fespace=&ifespace;
-
-      // set the parameters
-      plap_epsilon=new ConstantCoefficient(regularizationp);
-      plap_power=new ConstantCoefficient(powerp);
-      if (load==nullptr)
-      {
-         plap_input=new ConstantCoefficient(1.0);
-         input_ownership=true;
-      }
-      else
-      {
-         plap_input=load;
-         input_ownership=false;
-      }
-
-      nlform=nullptr;
-      nsolver=nullptr;
-      gmres=nullptr;
-      prec=nullptr;
-
-      // set the default integrator
-      integ=IntegratorType::HandCodedIntegrator;
-   }
-
-   ~ParNLSolverPLaplacian()
-   {
-      delete nlform;
-      delete nsolver;
-      delete prec;
-      delete gmres;
-      if (input_ownership) { delete plap_input;}
-      delete plap_epsilon;
-      delete plap_power;
-   }
-
-   /// Set the integrator.
-   /// 0 - hand coded, 1 - AD based (compute only Hessian by AD),
-   /// 2 - AD based (compute residual and Hessian by AD)
-   void SetIntegrator(IntegratorType intr)
-   {
-      integ=intr;
-   }
-
-   // set relative tolerance for the Newton solver
-   void SetNRRTol(real_t rtol)
-   {
-      newton_rtol=rtol;
-   }
-
-   // set absolute tolerance for the Newton solver
-   void SetNRATol(real_t atol)
-   {
-      newton_atol=atol;
-   }
-
-   // set max iterations for the NR solver
-   void SetMaxNRIter(int miter)
-   {
-      newton_iter=miter;
-   }
-
-   void SetLSRTol(real_t rtol)
-   {
-      linear_rtol=rtol;
-   }
-
-   void SetLSATol(real_t atol)
-   {
-      linear_atol=atol;
-   }
-
-   // set max iterations for the linear solver
-   void SetMaxLSIter(int miter)
-   {
-      linear_iter=miter;
-   }
-
-   // set the print level
-   void SetPrintLevel(int plev)
-   {
-      print_level=plev;
-   }
-
-   /// The state vector is used as initial condition for the NR solver. On
-   /// return the statev holds the solution to the problem.
-   void Solve(Vector& statev)
-   {
-      if (nlform==nullptr)
-      {
-         AllocSolvers();
-      }
-      Vector b; // RHS is zero
-      nsolver->Mult(b, statev);
-   }
-
-   /// Compute the energy
-   real_t GetEnergy(Vector& statev)
-   {
-      if (nlform==nullptr)
-      {
-         // allocate the solvers
-         AllocSolvers();
-      }
-      return nlform->GetEnergy(statev);
-   }
-
-private:
-   void AllocSolvers()
-   {
-      if (nlform!=nullptr) { delete nlform;}
-      if (nsolver!=nullptr) { delete nsolver;}
-      if (gmres!=nullptr) { delete gmres;}
-      if (prec!=nullptr) { delete prec;}
-
-      // Define the essential boundary attributes
-      Array<int> ess_bdr(mesh->bdr_attributes.Max());
-      ess_bdr = 1;
-
-      nlform = new ParNonlinearForm(fespace);
-      if (integ==IntegratorType::HandCodedIntegrator)
-      {
-         nlform->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,
-                                                  *plap_input));
-      }
-      else if (integ==IntegratorType::ADJacobianIntegrator)
-      {
-         // The template integrator is based on automatic differentiation. For
-         // ADJacobianIntegrator the residual (vector function) at an
-         // integration point is implemented as a functor by MyResidualFunctor.
-         // The vector function has a return size of four(4), four state
-         // arguments, and three(3) parameters. MyResidualFunctor is a template
-         // argument to the actual template class performing the differentiation
-         // - in this case, QVectorFuncAutoDiff. The derivatives are used in the
-         // integration loop in the integrator pLaplaceAD.
-         nlform->AddDomainIntegrator(new
-                                     pLaplaceAD<mfem::QVectorFuncAutoDiff<MyResidualFunctor,4,4,3>>(*plap_power,
-                                           *plap_epsilon,*plap_input));
-      }
-      else if (integ==IntegratorType::ADHessianIntegrator)
-      {
-         // The main difference from the previous case is that the user has to
-         // implement only a functional evaluation at an integration point. The
-         // implementation is in MyEnergyFunctor, which takes four state
-         // arguments and three parameters. The residual vector is the first
-         // derivative of the energy/functional with respect to the state
-         // variables, and the Hessian is the second derivative. Automatic
-         // differentiation is used for evaluating both of them.
-         nlform->AddDomainIntegrator(new
-                                     pLaplaceAD<mfem::QFunctionAutoDiff<MyEnergyFunctor,4,3>>(*plap_power,
-                                           *plap_epsilon,*plap_input));
-      }
-
-      nlform->SetEssentialBC(ess_bdr);
-
-      prec = new HypreBoomerAMG();
-      prec->SetPrintLevel(print_level);
-
-      gmres = new GMRESSolver(lcomm);
-      gmres->SetAbsTol(linear_atol);
-      gmres->SetRelTol(linear_rtol);
-      gmres->SetMaxIter(linear_iter);
-      gmres->SetPrintLevel(print_level);
-      gmres->SetPreconditioner(*prec);
-
-      nsolver = new NewtonSolver(lcomm);
-
-      nsolver->iterative_mode = true;
-      nsolver->SetSolver(*gmres);
-      nsolver->SetOperator(*nlform);
-      nsolver->SetPrintLevel(print_level);
-      nsolver->SetRelTol(newton_rtol);
-      nsolver->SetAbsTol(newton_atol);
-      nsolver->SetMaxIter(newton_iter);
-   }
-
-   real_t newton_rtol;
-   real_t newton_atol;
-   int newton_iter;
-
-   real_t linear_rtol;
-   real_t linear_atol;
-   int linear_iter;
-
-   int print_level;
-
-   // power of the p-laplacian
-   Coefficient* plap_power;
-   // regularization parameter
-   Coefficient* plap_epsilon;
-   // load(input) parameter
-   Coefficient* plap_input;
-   // flag indicating the ownership of plap_input
-   bool input_ownership;
-
-   MPI_Comm lcomm;
-
-   ParMesh *mesh;
-   ParFiniteElementSpace *fespace;
-
-   ParNonlinearForm *nlform;
-
-   HypreBoomerAMG *prec;
-   GMRESSolver *gmres;
-   NewtonSolver *nsolver;
-   IntegratorType integ;
-
-};
-
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI and HYPRE.
-   Mpi::Init(argc, argv);
-   int myrank = Mpi::WorldRank();
-   Hypre::Init();
-   // Define Caliper ConfigManager
-#ifdef MFEM_USE_CALIPER
-   cali::ConfigManager mgr;
-#endif
-   // Caliper instrumentation
-   MFEM_PERF_FUNCTION;
-
-   // 2. Parse command-line options
-   const char *mesh_file = "../../data/beam-tet.mesh";
-   int ser_ref_levels = 3;
-   int par_ref_levels = 1;
-   int order = 1;
-   bool visualization = true;
-   real_t newton_rel_tol = 1e-4;
-   real_t newton_abs_tol = 1e-6;
-   int newton_iter = 10;
-   int print_level = 0;
-
-   real_t pp = 2.0; // p-Laplacian power
-
-   IntegratorType integrator = IntegratorType::ADHessianIntegrator;
-   int int_integrator = integrator;
-   // HandCodedIntegrator  = 0 - do not use AD (hand coded)
-   // ADJacobianIntegrator = 1 - use AD for Hessian only
-   // ADHessianIntegrator  = 2 - use AD for Residual and Hessian
-
-   const char* cali_config = "runtime-report";
-
-   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(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&newton_rel_tol,
-                  "-rel",
-                  "--relative-tolerance",
-                  "Relative tolerance for the Newton solve.");
-   args.AddOption(&newton_abs_tol,
-                  "-abs",
-                  "--absolute-tolerance",
-                  "Absolute tolerance for the Newton solve.");
-   args.AddOption(&newton_iter,
-                  "-it",
-                  "--newton-iterations",
-                  "Maximum iterations for the Newton solve.");
-   args.AddOption(&pp,
-                  "-pp",
-                  "--power-parameter",
-                  "Power parameter (>=2.0) for the p-Laplacian.");
-   args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
-   args.AddOption(&int_integrator,
-                  "-int",
-                  "--integrator",
-                  "Integrator 0: standard; 1: AD for Hessian; 2: AD for residual and Hessian");
-   args.AddOption(&cali_config, "-p", "--caliper",
-                  "Caliper configuration string.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myrank == 0)
-      {
-         args.PrintUsage(std::cout);
-      }
-      return 1;
-   }
-   if (myrank == 0)
-   {
-      args.PrintOptions(std::cout);
-   }
-   integrator = static_cast<IntegratorType>(int_integrator);
-
-   StopWatch *timer = new StopWatch();
-
-   // Caliper configuration
-#ifdef MFEM_USE_CALIPER
-   mgr.add(cali_config);
-   mgr.start();
-#endif
-   // 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 load for the p-Laplacian
-   ConstantCoefficient load(1.00);
-
-   // 7. Define the finite element spaces for the solution
-   H1_FECollection fec(order, dim);
-   ParFiniteElementSpace fespace(pmesh, &fec, 1, Ordering::byVDIM);
-   HYPRE_Int glob_size = fespace.GlobalTrueVSize();
-   if (myrank == 0)
-   {
-      std::cout << "Number of finite element unknowns: " << glob_size
-                << std::endl;
-   }
-
-   // 8. Define the solution vector x as a parallel finite element grid function
-   //    corresponding to fespace. Initialize x with initial guess of zero,
-   //    which satisfies the boundary conditions.
-   ParGridFunction x(&fespace);
-   x = 0.0;
-   HypreParVector *sv = x.GetTrueDofs();
-
-   // 9. Define ParaView DataCollection
-   ParaViewDataCollection *dacol = new ParaViewDataCollection("Example",
-                                                              pmesh);
-   dacol->SetLevelsOfDetail(order);
-   dacol->RegisterField("sol", &x);
-
-   // 10. Define the NR solver
-   ParNLSolverPLaplacian* nr;
-
-   // 11. Start with linear diffusion - solvable for any initial guess
-   nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, 2.0, &load);
-   nr->SetIntegrator(integrator);
-   nr->SetMaxNRIter(newton_iter);
-   nr->SetNRATol(newton_abs_tol);
-   nr->SetNRRTol(newton_rel_tol);
-   nr->SetPrintLevel(print_level);
-   timer->Clear();
-   timer->Start();
-   nr->Solve(*sv);
-   timer->Stop();
-   if (myrank==0)
-   {
-      std::cout << "[pp=2] The solution time is: " << timer->RealTime()
-                << std::endl;
-   }
-   // Compute the energy
-   real_t energy = nr->GetEnergy(*sv);
-   if (myrank==0)
-   {
-      std::cout << "[pp=2] The total energy of the system is E=" << energy
-                << std::endl;
-   }
-   delete nr;
-   x.SetFromTrueDofs(*sv);
-   dacol->SetTime(2.0);
-   dacol->SetCycle(2);
-   dacol->Save();
-
-   // 12. Continue with powers higher than 2
-   for (int i = 3; i < pp; i++)
-   {
-      nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, (real_t)i, &load);
-      nr->SetIntegrator(integrator);
-      nr->SetMaxNRIter(newton_iter);
-      nr->SetNRATol(newton_abs_tol);
-      nr->SetNRRTol(newton_rel_tol);
-      nr->SetPrintLevel(print_level);
-      timer->Clear();
-      timer->Start();
-      nr->Solve(*sv);
-      timer->Stop();
-      if (myrank==0)
-      {
-         std::cout << "[pp="<<i<<"] The solution time is: " << timer->RealTime()
-                   << std::endl;
-      }
-      // Compute the energy
-      energy = nr->GetEnergy(*sv);
-      if (myrank==0)
-      {
-         std::cout << "[pp="<<i<<"] The total energy of the system is E=" << energy
-                   << std::endl;
-      }
-      delete nr;
-      x.SetFromTrueDofs(*sv);
-      dacol->SetTime((real_t)i);
-      dacol->SetCycle(i);
-      dacol->Save();
-   }
-
-   // 13. Continue with the final power
-   if (std::abs(pp - 2.0) > std::numeric_limits<real_t>::epsilon())
-   {
-      nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, pp, &load);
-      nr->SetIntegrator(integrator);
-      nr->SetMaxNRIter(newton_iter);
-      nr->SetNRATol(newton_abs_tol);
-      nr->SetNRRTol(newton_rel_tol);
-      nr->SetPrintLevel(print_level);
-      timer->Clear();
-      timer->Start();
-      nr->Solve(*sv);
-      timer->Stop();
-      if (myrank==0)
-      {
-         std::cout << "[pp="<<pp<<"] The solution time is: " << timer->RealTime()
-                   << std::endl;
-      }
-      // Compute the energy
-      energy = nr->GetEnergy(*sv);
-      if (myrank==0)
-      {
-         std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
-                   << std::endl;
-      }
-      delete nr;
-      x.SetFromTrueDofs(*sv);
-      dacol->SetTime(pp);
-      if (pp < 2.0)
-      {
-         dacol->SetCycle(static_cast<int>(std::floor(pp)));
-      }
-      else
-      {
-         dacol->SetCycle(static_cast<int>(std::ceil(pp)));
-      }
-      dacol->Save();
-   }
-
-   // 14. Free the used memory
-   delete dacol;
-   delete sv;
-   delete pmesh;
-   delete timer;
-
-   // Flush output before MPI_finalize
-#ifdef MFEM_USE_CALIPER
-   mgr.flush();
-#endif
-
-   return 0;
-}

miniapps/autodiff/pg.cpp

@@ -0,0 +1,78 @@
+#include "pg.hpp"
+namespace mfem
+{
+PGStepSizeRule::PGStepSizeRule(int rule_type,
+                               real_t alpha0, real_t max_alpha,
+                               real_t ratio, real_t ratio2)
+   : rule_type(static_cast<RuleType>(rule_type))
+   , max_alpha(max_alpha), alpha0(alpha0), ratio(ratio), ratio2(ratio2)
+{
+   MFEM_VERIFY(rule_type < RuleType::INVALID,
+               "PGStepSizeRule: Invalid rule type");
+   MFEM_VERIFY(alpha0 > 0, "PGStepSizeRule: alpha0 must be positive");
+   MFEM_VERIFY(max_alpha >= alpha0,
+               "PGStepSizeRule: max_alpha must be greater than or equal to alpha0");
+   if (rule_type == RuleType::CONSTANT)
+   {
+   }
+   else if (rule_type == RuleType::POLY)
+   {
+      MFEM_VERIFY(ratio > 0, "PGStepSizeRule: ratio must be positive for POLY rule");
+   }
+   else if (rule_type == RuleType::EXP)
+   {
+      MFEM_VERIFY(ratio > 1,
+                  "PGStepSizeRule: ratio must be greater than 1 for EXP rule");
+   }
+   else if (rule_type == RuleType::DOUBLE_EXP)
+   {
+      MFEM_VERIFY(ratio > 1 && ratio2 > 1,
+                  "PGStepSizeRule: ratio and ratio2 must be greater than 1 for DOUBLE_EXP rule");
+   }
+}
+
+real_t PGStepSizeRule::Get(int iter) const
+{
+   real_t alpha = alpha0;
+   switch (rule_type)
+   {
+      case RuleType::CONSTANT:
+         break;
+      case RuleType::POLY:
+         alpha *= std::pow(iter+1, ratio);
+         break;
+      case RuleType::EXP:
+         alpha *= std::pow(ratio, iter);
+         break;
+      case RuleType::DOUBLE_EXP:
+         alpha *= std::pow(ratio, std::pow(ratio2, iter));
+         break;
+      default:
+         break;
+   }
+   return std::min(alpha, max_alpha);
+}
+
+const GridFunction& ADPGFunctional::GetPrevLatent(int i) const
+{
+   Evaluator::param_t param = evaluator.Get(i);
+   const GridFunction* gf = std::visit([&](auto arg)
+   {
+      using T = std::decay_t<decltype(arg)>;
+      if constexpr (std::is_same_v<T, const GridFunction*>)
+      {
+         return (const GridFunction*)arg;
+      }
+      else
+      {
+         MFEM_ABORT("Parameter at index " << i
+                    << " is not a GridFunction or ParGridFunction");
+         return (const GridFunction*)nullptr;
+      }
+   }, param);
+   MFEM_VERIFY(gf != nullptr,
+               "ADPGFunctional: GetPrevLatent(" << i << ") is null");
+   return *gf;
+}
+
+}

miniapps/autodiff/pg.hpp

@@ -0,0 +1,472 @@
+#pragma once
+#include "mfem.hpp"
+#include "ad_native.hpp"
+#include "tools.hpp"
+
+namespace mfem
+{
+// PGStepSizeRule defines the step size rule for the Proximal Galerkin (PG) method.
+// See RuleType for the available rules
+struct PGStepSizeRule
+{
+   enum RuleType
+   {
+      CONSTANT, // alpha0
+      POLY, // alpha0 * (iter+1)^ratio
+      EXP, // alpha0 * ratio^iter
+      DOUBLE_EXP, // alpha0 * ratio^(ratio2^iter)
+      // ... add more rules as needed
+      INVALID // used to check for valid rule types
+   };
+   RuleType rule_type;
+
+   real_t max_alpha;
+   real_t alpha0; // initial step size
+   real_t ratio; // poly degree (POLY), exponential base (EXP, DOUBLE_EXP)
+   real_t ratio2; // nested exponential base (DOUBLE_EXP)
+
+   PGStepSizeRule(int rule_type,
+                  real_t alpha0 = 1.0, real_t max_alpha = 1e06,
+                  real_t ratio = -1.0, real_t ratio2 = -1.0);
+
+   /// Get the step size for the given iteration
+   real_t Get(int iter) const;
+};
+
+// Base struct for dual entropy functions
+class ADEntropy : public ADFunction
+{
+public:
+   ADEntropy(int n_input)
+      : ADFunction(n_input) { }
+   ADEntropy(int n_input, int capacity)
+      : ADFunction(n_input, capacity) { }
+};
+
+template <typename T>
+std::vector<T*> uniquevec2ptrvec(std::vector<std::unique_ptr<T>> &vec)
+{
+   std::vector<T*> ptrs(vec.size());
+   for (int i=0; i<vec.size(); i++)
+   {
+      ptrs[i] = vec[i].get();
+   }
+   return ptrs;
+}
+
+
+// Construct augmented energy for proximal Galerkin
+// psi =
+// L(u, psi) = f(u) + (1/alpha)(u*(psi-psi_k) - E^*(psi))
+// Equivalently, L(u, lambda) = f(u) + (u*lambda - E^*(alpha*lambda + psi_k))
+// so that
+// dL/du = df/du + (1/alpha)(psi-psi_k)
+// dL/dpsi = (1/alpha)(u - dE^*(psi))
+// When primal is not full vector, set primal_begin
+// The parameter should be [org_param, entropy_param, alpha, psi_k]
+class ADPGFunctional : public ADFunction
+{
+protected:
+   ADFunction &f;
+   std::vector<ADEntropy*> dual_entropy;
+   std::vector<int> primal_idx;
+   std::vector<int> dual_idx;
+   std::vector<int> entropy_size;
+   mutable const BlockVector *latent_k;
+   mutable Vector jac;
+   mutable DenseMatrix hess;
+   mutable real_t alpha;
+   std::unique_ptr<VectorCoefficient> owned_cf;
+   static int GetEntropySize(const std::vector<ADEntropy*> &dual_entropy)
+   {
+      int size = 0;
+      for (const auto &entropy : dual_entropy)
+      {
+         size += entropy->n_input;
+      }
+      return size;
+   }
+
+public:
+   ADPGFunctional(ADFunction &f, ADEntropy &dual_entropy, Evaluator::param_t alpha,
+                  int idx=0)
+      : ADFunction(f.n_input + dual_entropy.n_input, 1)
+      , f(f), dual_entropy{&dual_entropy}
+      , primal_idx(1)
+      , dual_idx(1)
+      , entropy_size(1)
+   {
+      evaluator.Add(alpha);
+      this->primal_idx[0] = idx;
+      entropy_size[0] = dual_entropy.n_input;
+      MFEM_VERIFY(f.n_input >= this->primal_idx[0] + entropy_size[0],
+                  "ADPGFunctional: f.n_input must not exceed "
+                  "primal_begin + dual_entropy.n_input:"
+                  << f.n_input << " >= " << n_input);
+      dual_idx[0] = f.n_input;
+   }
+   ADPGFunctional(ADFunction &f, ADEntropy &dual_entropy,
+                  Evaluator::param_t alpha,
+                  GridFunction &latent_k, int idx=0)
+      : ADPGFunctional(f, dual_entropy, alpha, idx)
+   {
+      evaluator.Add(&latent_k);
+   }
+   // Multiple entropies
+   ADPGFunctional(ADFunction &f, std::vector<ADEntropy*> dual_entropy_,
+                  std::vector<int> &primal_begin, Evaluator::param_t alpha)
+      : ADFunction(f.n_input + GetEntropySize(dual_entropy_), 1)
+      , f(f), dual_entropy(std::move(dual_entropy_))
+      , primal_idx(primal_begin)
+      , dual_idx(dual_entropy.size())
+      , entropy_size(dual_entropy.size())
+      , alpha(*evaluator.val.GetBlock(0).GetData())
+   {
+      evaluator.Add(alpha);
+      int dual_entropy_size = 0;
+      int max_primal_index = 0;
+      for (int i=0; i<dual_entropy.size(); i++)
+      {
+         dual_entropy_size += dual_entropy[i]->n_input;
+         max_primal_index = std::max(max_primal_index,
+                                     primal_begin[i] + dual_entropy[i]->n_input);
+      }
+      MFEM_VERIFY(f.n_input >= max_primal_index,
+                  "ADPGFunctional: f.n_input must be larger than "
+                  "primal_begin[i] + dual_entropy.n_input[i] for all i");
+   }
+
+   ADPGFunctional(ADFunction &f, std::vector<ADEntropy*> dual_entropy,
+                  std::vector<GridFunction*> latent_k_gf, std::vector<int> &primal_begin,
+                  Evaluator::param_t alpha)
+      : ADPGFunctional(f, std::move(dual_entropy), primal_begin, alpha)
+   {
+      MFEM_VERIFY(latent_k_gf.size() == this->dual_entropy.size(),
+                  "ADPGFunctional: latent_k must have the same size as dual_entropy: "
+                  << latent_k_gf.size() << " != " << dual_entropy.size());
+      MFEM_VERIFY(latent_k_gf.size() == primal_begin.size(),
+                  "ADPGFunctional: latent_k must have the same size as primal_begin"
+                  << latent_k_gf.size() << " != " << primal_begin.size());
+      for (int i=0; i<latent_k_gf.size(); i++)
+      {
+         MFEM_VERIFY(latent_k_gf[i] != nullptr,
+                     "ADPGFunctional: latent_k_gf[" << i << "] is null");
+         evaluator.Add(latent_k_gf[i]);
+      }
+   }
+   // Multiple entropies
+   ADPGFunctional(ADFunction &f,
+                  std::vector<std::unique_ptr<ADEntropy>> &dual_entropy,
+                  std::vector<int> &primal_begin, Evaluator::param_t alpha)
+      : ADPGFunctional(f, uniquevec2ptrvec(dual_entropy), primal_begin, alpha)
+   {}
+   ADPGFunctional(ADFunction &f,
+                  std::vector<std::unique_ptr<ADEntropy>> &dual_entropy,
+                  std::vector<std::unique_ptr<GridFunction>> &latent_k_gf,
+                  std::vector<int> primal_begin, Evaluator::param_t alpha)
+      : ADPGFunctional(f, uniquevec2ptrvec(dual_entropy),
+                       uniquevec2ptrvec(latent_k_gf), primal_begin, alpha)
+   {}
+
+   const GridFunction& GetPrevLatent(int i) const;
+
+   ADFunction &GetObjective() const
+   { return f; }
+
+   ADEntropy &GetEntropy() const
+   {
+      MFEM_VERIFY(dual_entropy.size() == 1,
+                  "ADPGFunctional: GetEntropy() can only be called when there is a single entropy");
+      return *dual_entropy[0];
+   }
+   const std::vector<ADEntropy*> &GetEntropies() const
+   { return dual_entropy; }
+
+   real_t GetAlpha() const { return alpha; }
+
+   void ProcessParameters(ElementTransformation &Tr,
+                          const IntegrationPoint &ip) const override
+   {
+      for (int i=0; i<dual_entropy.size(); i++)
+      {
+         dual_entropy[i]->ProcessParameters(Tr, ip);
+      }
+      f.ProcessParameters(Tr, ip);
+      latent_k = &evaluator.Eval(Tr, ip);
+      alpha = evaluator.val[0];
+   }
+
+   AD_IMPL(T, V, M, x_psi,
+   {
+      // variables
+      const V x(x_psi.GetData(), f.n_input);
+      V psi;
+
+      // evaluate mixed value
+      T cross_entropy = T();
+      T dual_entropy_sum = T();
+      for (int i=0; i<entropy_size.size(); i++)
+      {
+         psi.SetDataAndSize(x_psi.GetData() + dual_idx[i], entropy_size[i]);
+         const Vector &psi_k = latent_k->GetBlock(i+1);
+         for (int j=0; j<entropy_size[i]; j++)
+         {
+            cross_entropy += x[primal_idx[i] + j]*(psi[j] - psi_k[j]);
+         }
+         dual_entropy_sum += (*dual_entropy[i])(psi);
+      }
+      return f(x) + (cross_entropy - dual_entropy_sum)/alpha;
+   });
+};
+
+class ADLambdaPGFunctional : public ADPGFunctional
+{
+   using ADPGFunctional::ADPGFunctional;
+
+   AD_IMPL(T, V, M, x_lambda,
+   {
+      // variables
+      const V x(x_lambda.GetData(), f.n_input);
+      V lambda;
+      V latent;
+
+      // evaluate mixed value
+      T cross_entropy = T();
+      T dual_entropy_sum = T();
+      for (int i=0; i<entropy_size.size(); i++)
+      {
+         lambda.SetDataAndSize(x_lambda.GetData() + dual_idx[i], entropy_size[i]);
+         for (int j=0; j<entropy_size[i]; j++)
+         {
+            cross_entropy += x[primal_idx[i] + j]*lambda[j];
+         }
+         latent = latent_k->GetBlock(i+1);
+         latent.Add(alpha, lambda);
+         dual_entropy_sum += (*dual_entropy[i])(latent);
+      }
+      return f(x) + cross_entropy - dual_entropy_sum/alpha;
+   });
+};
+
+enum LatentType
+{
+   COEFFICIENT,
+   GF,
+   QF
+};
+
+
+// Dual entropy for (negative) Shannon entropy (xlogx - x) with half bound
+// when bound[1] = 1, [lower, inf[
+// when bound[1] = -1, ]-inf, upper]
+//
+// The resulting dual is (f(pm1*(x - shift)))^*
+// = f^*(pm1*x^*) + shift*pm1*x^*
+class ShannonEntropy : public ADEntropy
+{
+protected:
+   const real_t &bound;
+   int sign;
+public:
+   ShannonEntropy(Evaluator::param_t bound, int sign=1)
+      : ADEntropy(1, 1)
+      , bound(*evaluator.val.GetData())
+      , sign(sign)
+   {
+      evaluator.Add(bound);
+      MFEM_VERIFY(sign == 1 || sign == -1,
+                  "ShannonEntropy: sign must be 1 or -1");
+      MFEM_VERIFY(evaluator.val.GetBlock(0).Size() == 1,
+                  "ShannonEntropy: The provided bound has the wrong size. "
+                  "Expected 1, got " << evaluator.val.GetBlock(0).Size());
+   }
+   AD_IMPL(T, V, M, x, return sign*(exp(x[0]*sign)) + bound*x[0]; );
+};
+
+// Dual entropy for (negative) Fermi-Dirac with [lower, upper] bounds
+class FermiDiracEntropy : public ADEntropy
+{
+protected:
+   const real_t &upper_bound;
+   const real_t &lower_bound;
+   mutable real_t shift;
+   mutable real_t scale;
+public:
+   FermiDiracEntropy(Evaluator::param_t lower_bound,
+                     Evaluator::param_t upper_bound)
+      : ADEntropy(1, 2)
+      , upper_bound(*evaluator.val.GetData())
+      , lower_bound(*(evaluator.val.GetData()+1))
+   {
+      evaluator.Add(lower_bound);
+      evaluator.Add(upper_bound);
+      MFEM_VERIFY(evaluator.val.GetBlock(0).Size() == 1,
+                  "FermiDiracEntropy: The provided bound has the wrong size. "
+                  "Expected 1, got " << evaluator.val.GetBlock(0).Size());
+      MFEM_VERIFY(evaluator.val.GetBlock(1).Size() == 1,
+                  "FermiDiracEntropy: The provided bound has the wrong size. "
+                  "Expected 1, got " << evaluator.val.GetBlock(1).Size());
+   }
+   void ProcessParameters(const BlockVector &x) const override
+   {
+      shift = lower_bound;
+      scale = upper_bound - shift;
+   }
+   AD_IMPL(T, V, M, x,
+   {
+      T z = x[0]*scale;
+
+      // Use a numerically stable implementation of log(1+exp(z))
+      if (z > 0)
+      {
+         return z + log(1.0 + exp(-z)) + shift*x[0];
+      }
+      else
+      {
+         return log(1.0 + exp(z)) + shift*x[0];
+      }
+   });
+};
+// Dual entropy for (negative) Hellinger entropy with bound > 0
+class HellingerEntropy : public ADEntropy
+{
+   const real_t &scale;
+public:
+   HellingerEntropy(int dim, Evaluator::param_t bound)
+      : ADEntropy(dim, 1)
+      , scale(*evaluator.val.GetData())
+   {
+      evaluator.Add(bound);
+      MFEM_VERIFY(evaluator.val.GetBlock(0).Size() == 1,
+                  "HellingerEntropy: The provided bound has the wrong size. "
+                  "Expected 1, got " << evaluator.val.GetBlock(0).Size());
+   }
+   void ProcessParameters(const BlockVector &x) const override
+   {
+      MFEM_ASSERT(scale > 0, "HellingerEntropy: bound must be positive");
+   }
+   AD_IMPL(T, V, M, x, return sqrt(1 + (x*x)*(scale*scale)););
+};
+
+// Dual entropy for (negative) Simplex entropy with
+// x_i >= 0 sum_i x_i = bound
+// Also known as cateborical entropy or multinomial Shannon entropy
+class SimplexEntropy : public ADEntropy
+{
+   const real_t &scale;
+public:
+   SimplexEntropy(int n_input, Evaluator::param_t bound)
+      : ADEntropy(n_input, 1), scale(*evaluator.val.GetData())
+   {
+      evaluator.Add(bound);
+      MFEM_VERIFY(evaluator.val.GetBlock(0).Size() == 1,
+                  "SimplexEntropy: The provided bound has the wrong size. "
+                  "Expected 1, got " << evaluator.val.GetBlock(0).Size());
+   }
+
+   void ProcessParameters(const BlockVector &x) const override
+   {
+      MFEM_ASSERT(scale >= 0, "SimplexEntropy: bound must be non-negative");
+   }
+   AD_IMPL(T, V, M, x,
+   {
+      T maxval = x[0];
+      for (int i=1; i<x.Size(); i++) { maxval = max(maxval, x[i]); }
+
+      T sum_exp = T();
+      for (int i=0; i<x.Size(); i++)
+      {
+         sum_exp += exp(x[i]-maxval);
+      }
+      return scale*(maxval + log(sum_exp));
+   });
+};
+
+#ifdef MFEM_USE_PETSC
+class PetscOperatorWrapper : public Operator
+{
+protected:
+   MPI_Comm comm;
+   Operator &op;
+   Operator::Type mtype;
+   mutable std::unique_ptr<PetscParMatrix> petsc_matrix;
+public:
+   PetscOperatorWrapper(MPI_Comm comm, Operator &op,
+                        Operator::Type mtype = Operator::Type::PETSC_MATAIJ)
+      : Operator(op.Height(), op.Width()), comm(comm), op(op), mtype(mtype)
+   { }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      op.Mult(x, y);
+   }
+
+   Operator &GetGradient(const Vector &x) const override
+   {
+      auto &grad = op.GetGradient(x);
+      petsc_matrix = std::make_unique<PetscParMatrix>(comm, &grad, mtype);
+      return *petsc_matrix;
+   }
+};
+
+class NewtonLinearSolverMonitor : public IterativeSolverController
+{
+protected:
+   /// The last IterativeSolver to which this controller was attached.
+   const class IterativeSolver *iter_solver;
+#ifdef MFEM_USE_PETSC
+   PetscLinearSolver *petsc_solver;
+#endif
+   IterativeSolver *mfem_solver;
+
+   int numIterations=0;
+   int prefix=0;
+   bool is_root = true;
+   bool converged = false;
+
+public:
+#ifdef MFEM_USE_PETSC
+   NewtonLinearSolverMonitor(PetscLinearSolver &linear_solver)
+      : petsc_solver(&linear_solver)
+   {
+      is_root = Mpi::Root();
+   }
+#endif
+   NewtonLinearSolverMonitor(IterativeSolver &linear_solver)
+      : mfem_solver(&linear_solver)
+   {
+#ifdef MFEM_USE_MPI
+      is_root = Mpi::Root();
+#endif
+   }
+
+   void SetPrefix(size_t i) { prefix = i; }
+
+   virtual void Reset()
+   {
+      converged = false;
+      numIterations = 0;
+   }
+
+   /// Monitor the solution vector r
+   virtual void MonitorResidual(int it, real_t norm, const Vector &r,
+                                bool final)
+   {
+      if (final && is_root)
+      {
+         for (int i=0; i<prefix; i++) { out << " "; }
+         out << "Average Linear Solver Iterations: " << (numIterations /
+                                                         (it + 1.)) << std::endl;
+         numIterations = 0;
+         return;
+      }
+#ifdef MFEM_USE_PETSC
+      if (petsc_solver) { numIterations += petsc_solver->GetNumIterations(); }
+#endif
+      if (mfem_solver) { numIterations += mfem_solver->GetNumIterations(); }
+   }
+};
+#endif
+
+
+} // namespace mfem

miniapps/autodiff/seq_example.cpp

@@ -1,478 +0,0 @@
-// Copyright (c) 2010-2025, 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.
-//
-//                       MFEM AD Example - Serial Version
-//
-// Compile with: make seq_example
-//
-// Sample runs:  seq_example -m ../../data/beam-quad.mesh -pp 3.5
-//               seq_example -m ../../data/beam-tri.mesh  -pp 4.6
-//               seq_example -m ../../data/beam-hex.mesh
-//               seq_example -m ../../data/beam-tet.mesh
-//               seq_example -m ../../data/beam-wedge.mesh
-//
-// Description:  This examples solves a quasi-static nonlinear p-Laplacian
-//               problem with zero Dirichlet boundary conditions applied on all
-//               defined boundaries
-//
-//               The example demonstrates the use of nonlinear operators
-//               combined with automatic differentiation (AD). The integrators
-//               are defined in example.hpp. Selecting integrator = 0 will use
-//               the manually implemented integrator. Selecting integrator = 1
-//               or 2 will utilize one of the AD integrators.
-//
-//               We recommend viewing examples 1 and 19, before viewing this
-//               example.
-
-#include "example.hpp"
-
-using namespace mfem;
-
-enum IntegratorType
-{
-   HandCodedIntegrator  = 0,
-   ADJacobianIntegrator = 1,
-   ADHessianIntegrator  = 2
-};
-
-/// Non-linear solver for the p-Laplacian problem.
-class NLSolverPLaplacian
-{
-public:
-   /// Constructor Input: imesh - FE mesh, finite element space, power for the
-   /// p-Laplacian, external load (source, input), regularization parameter
-   NLSolverPLaplacian(Mesh& imesh, FiniteElementSpace& ifespace,
-                      real_t powerp=2,
-                      Coefficient* load=nullptr,
-                      real_t regularizationp=1e-7)
-   {
-      // default parameters for the Newton solver
-      newton_rtol = 1e-4;
-      newton_atol = 1e-8;
-      newton_iter = 10;
-
-      // linear solver
-      linear_rtol = 1e-7;
-      linear_atol = 1e-15;
-      linear_iter = 500;
-
-      print_level = 0;
-
-      // set the mesh
-      mesh=&imesh;
-
-      // set the fespace
-      fespace=&ifespace;
-
-      // set the parameters
-      plap_epsilon=new ConstantCoefficient(regularizationp);
-      plap_power=new ConstantCoefficient(powerp);
-      if (load==nullptr)
-      {
-         plap_input=new ConstantCoefficient(1.0);
-         input_ownership=true;
-      }
-      else
-      {
-         plap_input=load;
-         input_ownership=false;
-      }
-
-      // set the nonlinear form
-      nlform=nullptr;
-      lsolver=nullptr;
-      prec=nullptr;
-      nsolver=nullptr;
-
-      // set the default integrator
-      integ=IntegratorType::HandCodedIntegrator; // hand coded
-   }
-
-   ~NLSolverPLaplacian()
-   {
-      if (nlform!=nullptr) { delete nlform;}
-      if (nsolver!=nullptr) { delete nsolver;}
-      if (prec!=nullptr) { delete prec;}
-      if (lsolver!=nullptr) { delete lsolver;}
-      if (input_ownership) { delete plap_input;}
-      delete plap_epsilon;
-      delete plap_power;
-   }
-
-   /// Set the integrator.
-   /// 0 - hand coded, 1 - AD based (compute only Hessian by AD),
-   /// 2 - AD based (compute residual and Hessian by AD)
-   void SetIntegrator(IntegratorType intr)
-   {
-      integ=intr;
-   }
-
-
-   // set relative tolerance for the Newton solver
-   void SetNRRTol(real_t rtol)
-   {
-      newton_rtol=rtol;
-   }
-
-   // set absolute tolerance for the Newton solver
-   void SetNRATol(real_t atol)
-   {
-      newton_atol=atol;
-   }
-
-   // set max iterations for the NR solver
-   void SetMaxNRIter(int miter)
-   {
-      newton_iter=miter;
-   }
-
-   void SetLSRTol(real_t rtol)
-   {
-      linear_rtol=rtol;
-   }
-
-   void SetLSATol(real_t atol)
-   {
-      linear_atol=atol;
-   }
-
-   // set max iterations for the linear solver
-   void SetMaxLSIter(int miter)
-   {
-      linear_iter=miter;
-   }
-
-   // set the print level
-   void SetPrintLevel(int plev)
-   {
-      print_level=plev;
-   }
-
-   /// The state vector is used as initial condition for the NR solver. On
-   /// return the statev holds the solution to the problem.
-   void Solve(Vector& statev)
-   {
-      if (nlform==nullptr)
-      {
-         AllocSolvers();
-      }
-      Vector b; // RHS is zero
-      nsolver->Mult(b, statev);
-   }
-
-   /// Compute the energy
-   real_t GetEnergy(Vector& statev)
-   {
-      if (nlform==nullptr)
-      {
-         // allocate the solvers
-         AllocSolvers();
-      }
-      return nlform->GetEnergy(statev);
-   }
-
-
-private:
-
-   void AllocSolvers()
-   {
-      if (nlform!=nullptr) { delete nlform;}
-      if (nsolver!=nullptr) {delete nsolver;}
-      if (prec!=nullptr) {delete prec;}
-      if (lsolver!=nullptr) { delete lsolver;}
-
-      // Define the essential boundary attributes
-      Array<int> ess_bdr(mesh->bdr_attributes.Max());
-      ess_bdr = 1;
-
-      nlform = new NonlinearForm(fespace);
-
-      if (integ==IntegratorType::HandCodedIntegrator)
-      {
-         // standard hand coded integrator
-         nlform->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,
-                                                  *plap_input));
-      }
-      else if (integ==IntegratorType::ADJacobianIntegrator)
-      {
-         // The template integrator is based on automatic differentiation. For
-         // ADJacobianIntegrator the residual (vector function) at an
-         // integration point is implemented as a functor by MyVFunctor. The
-         // vector function has a return size of four(4), four state arguments,
-         // and three(3) parameters. MyVFunctor is a template argument to the
-         // actual template class performing the differentiation - in this case,
-         // QVectorFuncAutoDiff. The derivatives are used in the integration
-         // loop in the integrator pLaplaceAD.
-         nlform->AddDomainIntegrator(new
-                                     pLaplaceAD<mfem::QVectorFuncAutoDiff<MyResidualFunctor,4,4,3>>(*plap_power,
-                                           *plap_epsilon,*plap_input));
-      }
-      else // IntegratorType::ADHessianIntegrator
-      {
-         // The main difference from the previous case is that the user has to
-         // implement only a functional evaluation at an integration point. The
-         // implementation is in MyQFunctor, which takes four state arguments
-         // and three parameters. The residual vector is the first derivative of
-         // the energy/functional with respect to the state variables, and the
-         // Hessian is the second derivative. Automatic differentiation is used
-         // for evaluating both of them.
-         nlform->AddDomainIntegrator(new
-                                     pLaplaceAD<mfem::QFunctionAutoDiff<MyEnergyFunctor,4,3>>(*plap_power,
-                                           *plap_epsilon,*plap_input));
-      }
-
-      nlform->SetEssentialBC(ess_bdr);
-
-#ifdef MFEM_USE_SUITESPARSE
-      prec = new UMFPackSolver();
-#else
-      prec = new GSSmoother();
-#endif
-
-      // allocate the linear solver
-      lsolver=new CGSolver();
-      lsolver->SetRelTol(linear_rtol);
-      lsolver->SetAbsTol(linear_atol);
-      lsolver->SetMaxIter(linear_iter);
-      lsolver->SetPrintLevel(print_level);
-      lsolver->SetPreconditioner(*prec);
-
-      // allocate the NR solver
-      nsolver = new NewtonSolver();
-      nsolver->iterative_mode = true;
-      nsolver->SetSolver(*lsolver);
-      nsolver->SetOperator(*nlform);
-      nsolver->SetPrintLevel(print_level);
-      nsolver->SetRelTol(newton_rtol);
-      nsolver->SetAbsTol(newton_atol);
-      nsolver->SetMaxIter(newton_iter);
-   }
-
-   real_t newton_rtol;
-   real_t newton_atol;
-   int newton_iter;
-
-   real_t linear_rtol;
-   real_t linear_atol;
-   int linear_iter;
-
-   int print_level;
-
-   // reference to the mesh
-   Mesh* mesh;
-   // reference to the fespace
-   FiniteElementSpace *fespace;
-
-   // nonlinear form for the p-laplacian
-   NonlinearForm *nlform;
-   CGSolver *lsolver; // linear solver
-   Solver *prec; // preconditioner for the linear solver
-   NewtonSolver *nsolver; // NR solver
-   IntegratorType integ;
-
-   // power of the p-laplacian
-   Coefficient* plap_power;
-   // regularization parameter
-   Coefficient* plap_epsilon;
-   // load(input) parameter
-   Coefficient* plap_input;
-   // flag indicating the ownership of plap_input
-   bool input_ownership;
-};
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options
-   const char *mesh_file = "../../data/beam-tet.mesh";
-   int ser_ref_levels = 3;
-   int order = 1;
-   bool visualization = true;
-   real_t newton_rel_tol = 1e-4;
-   real_t newton_abs_tol = 1e-6;
-   int newton_iter = 10;
-   int print_level = 0;
-
-   real_t pp = 2.0; // p-Laplacian power
-
-   IntegratorType integrator = IntegratorType::ADHessianIntegrator;
-   int int_integrator = integrator;
-   // HandCodedIntegrator  = 0 - do not use AD (hand coded)
-   // ADJacobianIntegrator = 1 - use AD for Hessian only
-   // ADHessianIntegrator  = 2 - use AD for Residual and Hessian
-   StopWatch *timer = new StopWatch();
-   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(&order,
-                  "-o",
-                  "--order",
-                  "Order (degree) of the finite elements.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&newton_rel_tol,
-                  "-rel",
-                  "--relative-tolerance",
-                  "Relative tolerance for the Newton solve.");
-   args.AddOption(&newton_abs_tol,
-                  "-abs",
-                  "--absolute-tolerance",
-                  "Absolute tolerance for the Newton solve.");
-   args.AddOption(&newton_iter,
-                  "-it",
-                  "--newton-iterations",
-                  "Maximum iterations for the Newton solve.");
-   args.AddOption(&pp,
-                  "-pp",
-                  "--power-parameter",
-                  "Power parameter (>=2.0) for the p-Laplacian.");
-   args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
-   args.AddOption(&int_integrator,
-                  "-int",
-                  "--integrator",
-                  "Integrator 0: standard; 1: AD for Hessian; 2: AD for residual and Hessian");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(std::cout);
-      return 1;
-   }
-   args.PrintOptions(std::cout);
-   integrator = static_cast<IntegratorType>(int_integrator);
-
-   // 2. Read the (serial) mesh from the given mesh file.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 3. 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();
-   }
-
-   // 4. Define the load parameter for the p-Laplacian
-   ConstantCoefficient load(1.00);
-
-   // 5. Define the finite element spaces for the solution
-   H1_FECollection fec(order, dim);
-   FiniteElementSpace fespace(mesh, &fec, 1, Ordering::byVDIM);
-   int glob_size = fespace.GetTrueVSize();
-
-   std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
-
-   // 6. Define the solution grid function
-   GridFunction x(&fespace);
-   x = 0.0;
-
-   // 7. Define the solution true vector
-   Vector sv(fespace.GetTrueVSize());
-   sv = 0.0;
-
-   // 8. Define ParaView DataCollection
-   ParaViewDataCollection *dacol = new ParaViewDataCollection("Example", mesh);
-   dacol->SetLevelsOfDetail(order);
-   dacol->RegisterField("sol", &x);
-
-   // 9. Define the nonlinear p-Laplacian solver
-   NLSolverPLaplacian* nr;
-
-   // 10. Start with linear diffusion - solvable for any initial guess
-   nr=new NLSolverPLaplacian(*mesh, fespace, 2.0, &load);
-   nr->SetIntegrator(integrator);
-   nr->SetMaxNRIter(newton_iter);
-   nr->SetNRATol(newton_abs_tol);
-   nr->SetNRRTol(newton_rel_tol);
-   timer->Clear();
-   timer->Start();
-   nr->Solve(sv);
-   timer->Stop();
-   std::cout << "[pp=2] The solution time is: " << timer->RealTime()
-             << std::endl;
-   // Compute the energy
-   real_t energy = nr->GetEnergy(sv);
-   std::cout << "[pp=2] The total energy of the system is E=" << energy
-             << std::endl;
-   delete nr;
-   x.SetFromTrueDofs(sv);
-   dacol->SetTime(2.0);
-   dacol->SetCycle(2);
-   dacol->Save();
-
-
-   // 11. Continue with powers higher than 2
-   for (int i = 3; i < pp; i++)
-   {
-      nr=new NLSolverPLaplacian(*mesh, fespace, (real_t)i, &load);
-      nr->SetIntegrator(integrator);
-      nr->SetMaxNRIter(newton_iter);
-      nr->SetNRATol(newton_abs_tol);
-      nr->SetNRRTol(newton_rel_tol);
-      timer->Clear();
-      timer->Start();
-      nr->Solve(sv);
-      timer->Stop();
-      std::cout << "[pp=" << i
-                << "] The solution time is: " << timer->RealTime() << std::endl;
-      energy = nr->GetEnergy(sv);
-      std::cout << "[pp="<< i<<"] The total energy of the system is E=" << energy
-                << std::endl;
-      delete nr;
-      x.SetFromTrueDofs(sv);
-      dacol->SetTime(i);
-      dacol->SetCycle(i);
-      dacol->Save();
-   }
-
-   // 12. Continue with the final power
-   if (std::abs(pp - 2.0) > std::numeric_limits<real_t>::epsilon())
-   {
-      nr=new NLSolverPLaplacian(*mesh, fespace, pp, &load);
-      nr->SetIntegrator(integrator);
-      nr->SetMaxNRIter(newton_iter);
-      nr->SetNRATol(newton_abs_tol);
-      nr->SetNRRTol(newton_rel_tol);
-      timer->Clear();
-      timer->Start();
-      nr->Solve(sv);
-      timer->Stop();
-      std::cout << "[pp=" << pp
-                << "] The solution time is: " << timer->RealTime() << std::endl;
-      energy = nr->GetEnergy(sv);
-      std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
-                << std::endl;
-      delete nr;
-      x.SetFromTrueDofs(sv);
-      dacol->SetTime(pp);
-      if (pp < 2.0)
-      {
-         dacol->SetCycle(static_cast<int>(std::floor(pp)));
-      }
-      else
-      {
-         dacol->SetCycle(static_cast<int>(std::ceil(pp)));
-      }
-      dacol->Save();
-   }
-
-   // 13. Free the memory
-   delete dacol;
-   delete mesh;
-   delete timer;
-   return 0;
-}

miniapps/autodiff/seq_test.cpp

@@ -1,158 +0,0 @@
-// Copyright (c) 2010-2025, 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 "admfem.hpp"
-#include "mfem.hpp"
-
-template<typename TDataType, typename TParamVector, typename TStateVector
-         , int state_size, int param_size>
-class DiffusionFunctional
-{
-public:
-   TDataType operator() (TParamVector& vparam, TStateVector& uu)
-   {
-      MFEM_ASSERT(state_size==4,"ExampleFunctor state_size should be equal to 4!");
-      MFEM_ASSERT(param_size==2,"ExampleFunctor param_size should be equal to 2!");
-      auto kappa = vparam[0]; // diffusion coefficient
-      auto load = vparam[1];  // volumetric influx
-      TDataType rez = kappa*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])/2.0 - load*uu[3];
-      return rez;
-   }
-
-};
-
-template<typename TDataType, typename TParamVector, typename TStateVector,
-         int residual_size, int state_size, int param_size>
-class DiffusionResidual
-{
-public:
-   void operator ()(TParamVector& vparam, TStateVector& uu, TStateVector& rr)
-   {
-      MFEM_ASSERT(residual_size==4,
-                  "DiffusionResidual residual_size should be equal to 4!");
-      MFEM_ASSERT(state_size==4,"ExampleFunctor state_size should be equal to 4!");
-      MFEM_ASSERT(param_size==2,"ExampleFunctor param_size should be equal to 2!");
-      auto kappa = vparam[0]; // diffusion coefficient
-      auto load = vparam[1];  // volumetric influx
-
-      rr[0] = kappa * uu[0];
-      rr[1] = kappa * uu[1];
-      rr[2] = kappa * uu[2];
-      rr[3] = -load;
-   }
-};
-
-int main(int argc, char *argv[])
-{
-
-#ifdef MFEM_USE_ADFORWARD
-   std::cout<<"MFEM_USE_ADFORWARD == true"<<std::endl;
-#else
-   std::cout<<"MFEM_USE_ADFORWARD == false"<<std::endl;
-#endif
-
-#ifdef MFEM_USE_CALIPER
-   cali::ConfigManager mgr;
-#endif
-   // Caliper instrumentation
-   MFEM_PERF_FUNCTION;
-#ifdef MFEM_USE_CALIPER
-   const char* cali_config = "runtime-report";
-   mgr.add(cali_config);
-   mgr.start();
-#endif
-   mfem::Vector param(2);
-   param[0]=3.0; // diffusion coefficient
-   param[1]=2.0; // volumetric influx
-
-   mfem::Vector state(4);
-   state[0]=1.0; // grad_x
-   state[1]=2.0; // grad_y
-   state[2]=3.0; // grad_z
-   state[3]=4.0; // state value
-
-   mfem::QFunctionAutoDiff<DiffusionFunctional,4,2> adf;
-
-   mfem::Vector rr0(4);
-   mfem::DenseMatrix hh0(4,4);
-
-   mfem::Vector rr1(4);
-   mfem::DenseMatrix hh1(4,4);
-   MFEM_PERF_BEGIN("Grad");
-   adf.Grad(param,state,rr0);
-   MFEM_PERF_END("Grad");
-   MFEM_PERF_BEGIN("Hessian");
-   adf.Hessian(param, state, hh0);
-   MFEM_PERF_END("Hessian");
-   // dump out the results
-   std::cout<<"FunctionAutoDiff"<<std::endl;
-   std::cout<< adf.Eval(param,state)<<std::endl;
-   rr0.Print(std::cout);
-   hh0.Print(std::cout);
-
-   mfem::QVectorFuncAutoDiff<DiffusionResidual,4,4,2> rdf;
-   MFEM_PERF_BEGIN("Jacobian");
-   rdf.Jacobian(param, state, hh1);
-   MFEM_PERF_END("Jacobian");
-
-   std::cout<<"ResidualAutoDiff"<<std::endl;
-   hh1.Print(std::cout);
-
-   // using lambda expression
-   auto func = [](mfem::Vector& vparam,
-                  mfem::ad::ADVectorType& uu,
-                  mfem::ad::ADVectorType& vres)
-   {
-      // auto func = [](auto& vparam, auto& uu, auto& vres) { //c++14
-      auto kappa = vparam[0]; // diffusion coefficient
-      auto load = vparam[1];  // volumetric influx
-
-      vres[0] = kappa * uu[0];
-      vres[1] = kappa * uu[1];
-      vres[2] = kappa * uu[2];
-      vres[3] = -load;
-   };
-
-   mfem::VectorFuncAutoDiff<4,4,2> fdr(func);
-   MFEM_PERF_BEGIN("JacobianV");
-   fdr.Jacobian(param,state,
-                hh1); // computes the gradient of func and stores the result in hh1
-   MFEM_PERF_END("JacobianV");
-   std::cout<<"LambdaAutoDiff"<<std::endl;
-   hh1.Print(std::cout);
-
-
-   mfem::real_t kappa = param[0];
-   mfem::real_t load =  param[1];
-   // using lambda expression
-   auto func01 = [&kappa,&load](mfem::Vector& vparam,
-                                mfem::ad::ADVectorType& uu,
-                                mfem::ad::ADVectorType& vres)
-   {
-      // auto func = [](auto& vparam, auto& uu, auto& vres) { //c++14
-
-      vres[0] = kappa * uu[0];
-      vres[1] = kappa * uu[1];
-      vres[2] = kappa * uu[2];
-      vres[3] = -load;
-   };
-
-   mfem::VectorFuncAutoDiff<4,4,2> fdr01(func01);
-   MFEM_PERF_BEGIN("Jacobian1");
-   fdr01.Jacobian(param,state,hh1);
-   MFEM_PERF_END("Jacobian1");
-   std::cout<<"LambdaAutoDiff 01"<<std::endl;
-   hh1.Print(std::cout);
-
-#ifdef MFEM_USE_CALIPER
-   mgr.flush();
-#endif
-}

miniapps/autodiff/tadvector.hpp

@@ -16,7 +16,6 @@
 
 #include <cmath>
 #include <iostream>
-#include <limits>
 #if defined(_MSC_VER) && (_MSC_VER < 1800)
 #include <float.h>
 #define isfinite _finite
@@ -119,11 +118,7 @@ public:
        with SetData(). */
    TAutoDiffVector(dtype *_data, int _size)
    {
-      if (capacity > 0)
-      {
-         delete[] data;
-         capacity = 0;
-      }
+      capacity = 0;
       size = _size;
       data = _data;
    }
@@ -315,7 +310,7 @@ public:
    /// Dot product with a `dtype *` array.
    dtype operator*(const dtype *v) const
    {
-      dtype dot = 0.0;
+      dtype dot = {};
       for (int i = 0; i < size; i++)
       {
          dot += data[i] * v[i];
@@ -327,7 +322,7 @@ public:
    dtype operator*(const TAutoDiffVector<dtype> &v) const
    {
       MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
-      dtype dot = 0.0;
+      dtype dot = {};
       for (int i = 0; i < size; i++)
       {
          dot += data[i] * v[i];
@@ -338,7 +333,7 @@ public:
    dtype operator*(const Vector &v) const
    {
       MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
-      dtype dot = 0.0;
+      dtype dot = {};
       for (int i = 0; i < size; i++)
       {
          dot += data[i] * v[i];
@@ -586,7 +581,7 @@ public:
    }
 
    /// Destroys vector.
-   ~TAutoDiffVector() { delete[] data; }
+   ~TAutoDiffVector() { if (OwnsData()) { delete[] data; } }
 
    /// Prints vector to stream @a os with @a width entries per line.
    void Print(std::ostream &os = mfem::out, int width = 8) const
@@ -649,8 +644,8 @@ public:
          return abs(data[0]);
       } // end if 1 == size
 
-      dtype scale = 0.0;
-      dtype sum = 0.0;
+      dtype scale = {};
+      dtype sum = {};
 
       for (int i = 0; i < size; i++)
       {
@@ -674,7 +669,7 @@ public:
    /// Returns the l_infinity norm of the vector.
    dtype Normlinf() const
    {
-      dtype max = 0.0;
+      dtype max = {};
       for (int i = 0; i < size; i++)
       {
          max = max(abs(data[i]), max);
@@ -684,7 +679,7 @@ public:
    /// Returns the l_1 norm of the vector.
    dtype Norml1() const
    {
-      dtype sum = 0.0;
+      dtype sum = {};
       for (int i = 0; i < size; i++)
       {
          sum += abs(data[i]);

miniapps/autodiff/tools.hpp

@@ -0,0 +1,230 @@
+#pragma once
+#include "mfem.hpp"
+
+namespace mfem
+{
+class MappedGridFunctionCoefficient : public Coefficient
+{
+private:
+   GridFunction *gf;
+   std::function<real_t(const real_t)> map_func;
+public:
+   MappedGridFunctionCoefficient(GridFunction *gf_,
+                                 std::function<real_t(const real_t)> map_func_)
+      : gf(gf_), map_func(map_func_) {  }
+   virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      return map_func(gf->GetValue(T.ElementNo, T.GetIntPoint()));
+   }
+};
+class VectorGradientGridFunction : public MatrixCoefficient
+{
+private:
+   GridFunction &gf;
+public:
+   VectorGradientGridFunction(GridFunction &gf)
+      : MatrixCoefficient(gf.FESpace()->GetVDim(),
+                          gf.FESpace()->GetMesh()->SpaceDimension()), gf(gf)
+   {}
+
+   void Eval(DenseMatrix &grad, ElementTransformation &T,
+             const IntegrationPoint &ip) override
+   { gf.GetVectorGradient(T, grad); }
+};
+
+inline std::unique_ptr<GridFunction>
+NewGridFunction(FiniteElementSpace &fes)
+{
+#ifdef MFEM_USE_MPI
+   if (ParFiniteElementSpace *pfes =
+          dynamic_cast<ParFiniteElementSpace*>(&fes))
+   {
+      return std::make_unique<ParGridFunction>(pfes);
+   }
+#endif
+   return std::make_unique<GridFunction>(&fes);
+}
+
+inline std::unique_ptr<LinearForm>
+NewLinearForm(FiniteElementSpace &fes)
+{
+#ifdef MFEM_USE_MPI
+   if (ParFiniteElementSpace *pfes =
+          dynamic_cast<ParFiniteElementSpace*>(&fes))
+   {
+      return std::make_unique<ParLinearForm>(pfes);
+   }
+#endif
+   return std::make_unique<LinearForm>(&fes);
+}
+
+inline std::unique_ptr<BilinearForm>
+NewBilinearForm(FiniteElementSpace &fes)
+{
+#ifdef MFEM_USE_MPI
+   if (ParFiniteElementSpace *pfes =
+          dynamic_cast<ParFiniteElementSpace*>(&fes))
+   {
+      return std::make_unique<ParBilinearForm>(pfes);
+   }
+#endif
+   return std::make_unique<BilinearForm>(&fes);
+}
+
+inline std::unique_ptr<MixedBilinearForm>
+NewMixedBilinearForm(FiniteElementSpace &trial_fes,
+                     FiniteElementSpace &test_fes)
+{
+#ifdef MFEM_USE_MPI
+   if (ParFiniteElementSpace *trial_pfes =
+          dynamic_cast<ParFiniteElementSpace*>(&trial_fes))
+   {
+      ParFiniteElementSpace *test_pfes = dynamic_cast<ParFiniteElementSpace*>
+                                         (&test_fes);
+      MFEM_VERIFY(test_pfes != nullptr,
+                  "NewMixedBilinearForm: Trial is parallel, but test is not.");
+      return std::make_unique<ParMixedBilinearForm>(trial_pfes, test_pfes);
+   }
+   MFEM_VERIFY(dynamic_cast<ParFiniteElementSpace*>(&test_fes) == nullptr,
+               "NewMixedBilinearForm: Trial is not parallel, but test is.");
+#endif
+   return std::make_unique<MixedBilinearForm>(&trial_fes, &test_fes);
+}
+inline std::unique_ptr<NonlinearForm>
+NewNonlinearForm(FiniteElementSpace &fes)
+{
+#ifdef MFEM_USE_MPI
+   if (ParFiniteElementSpace *pfes =
+          dynamic_cast<ParFiniteElementSpace*>(&fes))
+   {
+      return std::make_unique<ParNonlinearForm>(pfes);
+   }
+#endif
+   return std::make_unique<NonlinearForm>(&fes);
+}
+inline std::unique_ptr<BlockNonlinearForm>
+NewBlockNonlinearForm(Array<FiniteElementSpace*> &fes)
+{
+#ifdef MFEM_USE_MPI
+   int numParallel = 0;
+
+   Array<ParFiniteElementSpace*> pfes;
+   for (auto *space : fes)
+   {
+      pfes.Append(dynamic_cast<ParFiniteElementSpace*>(space));
+      numParallel += pfes.Last() != nullptr;
+   }
+   MFEM_VERIFY(numParallel == 0 || numParallel == fes.Size(),
+               "NewBlockNonlinearForm: either all or none of the spaces must be parallel");
+   if (numParallel == fes.Size())
+   {
+      return std::make_unique<ParBlockNonlinearForm>(pfes);
+   }
+#endif
+   return std::make_unique<BlockNonlinearForm>(fes);
+}
+
+#ifdef MFEM_USE_PETSC
+// Monolithic direct solver for block system
+class MUMPSMonoSolver : public MUMPSSolver
+{
+private:
+   std::unique_ptr<HypreParMatrix> mono;
+public:
+   MUMPSMonoSolver(MPI_Comm comm) : MUMPSSolver(comm) {}
+
+   void SetOperator(const Operator &op)
+   {
+      const BlockOperator *bop = dynamic_cast<const BlockOperator*>(&op);
+      MFEM_VERIFY(bop != nullptr, "Not a BlockOperator");
+      Array2D<const HypreParMatrix*> blocks(bop->NumRowBlocks(), bop->NumColBlocks());
+      for (int j=0; j<bop->NumColBlocks(); j++)
+      {
+         for (int i=0; i<bop->NumRowBlocks(); i++)
+         {
+            if (bop->IsZeroBlock(i,j)) { continue; }
+            const HypreParMatrix *m =
+               dynamic_cast<const HypreParMatrix*>(&bop->GetBlock(i,j));
+            MFEM_VERIFY(m != nullptr, "Not a HypreParMatrix");
+            blocks(i,j) = m;
+         }
+      }
+      mono.reset(HypreParMatrixFromBlocks(blocks));
+      MUMPSSolver::SetOperator(*mono);
+   };
+};
+#endif
+
+inline std::tuple<std::unique_ptr<FiniteElementSpace>, std::unique_ptr<L2_FECollection>>
+      QSpaceToFESpace(QuadratureSpace &qs)
+{
+   Mesh *mesh = qs.GetMesh();
+   const int dim = mesh->Dimension();
+   Geometry::Type geom = mesh->GetTypicalElementGeometry();
+   MFEM_VERIFY(geom != Geometry::TRIANGLE &&
+               geom != Geometry::TETRAHEDRON &&
+               geom != Geometry::PRISM &&
+               geom != Geometry::PYRAMID,
+               "QSpaceToFESpace: only support tensor product elements");
+   std::unique_ptr<L2_FECollection> fec
+      = std::make_unique<L2_FECollection> (qs.GetOrder()/2, dim);
+
+   std::unique_ptr<FiniteElementSpace> fes;
+#ifdef MFEM_USE_MPI
+   ParMesh *pmesh = dynamic_cast<ParMesh*>(qs.GetMesh());
+   if (pmesh) { fes = std::make_unique<ParFiniteElementSpace>(pmesh, fec.get()); }
+#endif
+   if (!fes) { fes = std::make_unique<FiniteElementSpace>(mesh, fec.get()); }
+   return std::make_tuple(std::move(fes), std::move(fec));
+}
+
+inline Array<int> GetOffsets(const Array<FiniteElementSpace*> &fespaces)
+{
+   Array<int> offsets(fespaces.Size() + 1);
+   offsets[0] = 0;
+   for (int i=0; i<fespaces.Size(); i++)
+   {
+      offsets[i+1] = offsets[i] + fespaces[i]->GetVSize();
+   }
+   return std::move(offsets);
+}
+inline Array<int> GetTrueOffsets(const Array<FiniteElementSpace*> &fespaces)
+{
+   Array<int> offsets(fespaces.Size() + 1);
+   offsets[0] = 0;
+   for (int i=0; i<fespaces.Size(); i++)
+   {
+      offsets[i+1] = offsets[i] + fespaces[i]->GetTrueVSize();
+   }
+   return std::move(offsets);
+}
+
+class VectorNormCoefficient : public Coefficient
+{
+private:
+   VectorCoefficient &vc;
+   Vector v;
+public:
+   VectorNormCoefficient(VectorCoefficient &vc): vc(vc), v(vc.GetVDim()) {}
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      vc.Eval(v, T, ip);
+      return std::sqrt(v*v);
+   }
+};
+
+class BooleanCoefficient : public Coefficient
+{
+private:
+   Coefficient &cf;
+   std::function<bool(real_t)> func;
+public:
+   BooleanCoefficient(Coefficient &cf, std::function<bool(real_t)> func)
+      : cf(cf), func(func) {}
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      return func(cf.Eval(T, ip));
+   }
+};
+
+};