size error fix

fem/CMakeLists.txt

@@ -57,6 +57,7 @@ set(SRCS
   integ/nonlininteg_vecconvection_pa.cpp
   integ/nonlininteg_vecconvection_mf.cpp
   coefficient.cpp
+  autodiff.cpp
   complex_fem.cpp
   convergence.cpp
   datacollection.cpp
@@ -165,6 +166,7 @@ set(HDRS
   integ/bilininteg_hcurlhdiv_kernels.hpp
   integ/bilininteg_mass_kernels.hpp
   coefficient.hpp
+  autodiff.hpp
   complex_fem.hpp
   convergence.hpp
   datacollection.hpp

fem/autodiff.cpp

@@ -0,0 +1,335 @@
+
+// Copyright (c) 2010-2024, 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 "autodiff.hpp"
+
+namespace mfem
+{
+
+void ADVectorFunc::Jacobian(const Vector &vstate,DenseMatrix &jac)
+{
+   MFEM_ASSERT(vstate.Size() == jac.Width(),
+               "state and jacobian size do not match");
+
+   jac = 0.0;
+#ifdef MFEM_USE_CODIPACK
+   int vector_size = jac.Height();
+   int state_size = jac.Width();
+   ADVector ad_state(state_size);
+   ADVector ad_result(vector_size);
+#ifdef MFEM_USE_ADFORWARD
+   for (int i=0; i<state_size; i++)
+   {
+      ad_state[i].setValue(vstate[i]);
+      ad_state[i].setGradient(0.0);
+   }
+   for (int ii=0; ii<state_size; ii++)
+   {
+      ad_state[ii].setGradient(1.0);
+      F(param,ad_state,ad_result);
+      for (int jj=0; jj<vector_size; jj++)
+      {
+         jac(jj,ii)=ad_result[jj].getGradient();
+      }
+      ad_state[ii].setGradient(0.0);
+   }
+#else // use reverse mode
+   for (int i=0; i<state_size; i++)
+   {
+      ad_state[i]=vstate[i];
+   }
+
+   ADFloat::Tape& tape = ADFloat::getTape();
+   typename ADFloat::Tape::Position pos = tape.getPosition();
+
+   tape.setActive();
+   for (int ii=0; ii<state_size; ii++) { tape.registerInput(ad_state[ii]); }
+   F(param,ad_state,ad_result);
+   for (int ii=0; ii<vector_size; ii++) { tape.registerOutput(ad_result[ii]); }
+   tape.setPassive();
+
+   for (int jj=0; jj<vector_size; jj++)
+   {
+      ad_result[jj].setGradient(1.0);
+      tape.evaluate();
+      for (int ii=0; ii<state_size; ii++)
+      {
+         jac(jj,ii)=ad_state[ii].getGradient();
+      }
+      tape.clearAdjoints();
+      ad_result[jj].setGradient(0.0);
+   }
+   tape.resetTo(pos);
+#endif
+#else
+   MFEM_WARNING("not compiled with CoDiPack -- will return zero gradient");
+#endif
+}
+
+void ADVectorFunc::Curl(const Vector &vstate, Vector &curl)
+{
+   MFEM_ASSERT(vstate.Size() == 3, "state size expected to be 3");
+   MFEM_ASSERT(curl.Size() == 3, "solution size expected to be 3");
+   DenseMatrix jac(vstate.Size()); // Assume the jacobian is square
+   Jacobian(vstate, jac);
+
+   curl[0] = jac(2,1) - jac(1,2);
+   curl[1] = jac(0,2) - jac(2,0);
+   curl[2] = jac(1,0) - jac(0,1);
+}
+
+real_t ADVectorFunc::Divergence(const Vector &vstate)
+{
+   DenseMatrix jac(vstate.Size()); // Assume the jacobian is square
+   Jacobian(vstate, jac);
+   real_t div = 0.0;
+   for (int ii = 0; ii<vstate.Size(); ii++)
+   {
+      div += jac(ii,ii);
+   }
+   return div;
+}
+
+void ADVectorFunc::Gradient(const Vector &vstate, Vector &grad)
+{
+   MFEM_ASSERT( grad.Size() == vstate.Size(), "grad and state size not equal");
+
+   DenseMatrix jac(1, vstate.Size());
+   Jacobian(vstate, jac);
+
+   jac.GetRow(0, grad);
+}
+
+void ADVectorFunc::Solution(const Vector &vstate, Vector &sol)
+{
+   ADVector ad_state(vstate.Size());
+   ADVector ad_result(sol.Size());
+   for (int i=0; i<vstate.Size(); i++)
+   {
+#ifdef MFEM_USE_CODIPACK
+      ad_state[i].setValue(vstate[i]);
+      ad_state[i].setGradient(0.0);
+#else
+      ad_state[i] = vstate[i];
+#endif
+   }
+
+   F(param,ad_state,ad_result);
+   for (int jj=0; jj<sol.Size(); jj++)
+   {
+#ifdef MFEM_USE_CODIPACK
+      sol[jj] = ad_result[jj].getValue();
+#else
+      sol[jj] = ad_result[jj];
+#endif
+   }
+}
+
+real_t ADVectorFunc::ScalarSolution(const Vector &vstate)
+{
+   Vector sol(1);
+   Solution(vstate, sol);
+   return sol[0];
+}
+
+void ADVectorTDFunc::Jacobian(const Vector &vstate, const real_t time,
+                              DenseMatrix &jac)
+{
+   MFEM_ASSERT(vstate.Size() == jac.Width(),
+               "state and jacobian size do not match");
+
+   jac = 0.0;
+#ifdef MFEM_USE_CODIPACK
+   int vector_size = jac.Height();
+   int state_size = jac.Width();
+   ADVector ad_state(state_size);
+   ADVector ad_result(vector_size);
+   ADFloat ad_t(time);
+#ifdef MFEM_USE_ADFORWARD
+   ad_t.setGradient(0.0);
+   for (int i=0; i<state_size; i++)
+   {
+      ad_state[i].setValue(vstate[i]);
+      ad_state[i].setGradient(0.0);
+   }
+
+   for (int ii=0; ii<state_size; ii++)
+   {
+      ad_state[ii].setGradient(1.0);
+      F(param,ad_state, ad_t, ad_result);
+      for (int jj=0; jj<vector_size; jj++)
+      {
+         jac(jj,ii)=ad_result[jj].getGradient();
+      }
+      ad_state[ii].setGradient(0.0);
+   }
+#else // use reverse mode
+   for (int i=0; i<state_size; i++)
+   {
+      ad_state[i]=vstate[i];
+   }
+
+   ADFloat::Tape& tape = ADFloat::getTape();
+   typename ADFloat::Tape::Position pos = tape.getPosition();
+
+   tape.setActive();
+   for (int ii=0; ii<state_size; ii++) { tape.registerInput(ad_state[ii]); }
+   F(param,ad_state,ad_t,ad_result);
+   for (int ii=0; ii<vector_size; ii++) { tape.registerOutput(ad_result[ii]); }
+   tape.setPassive();
+
+   for (int jj=0; jj<vector_size; jj++)
+   {
+      ad_result[jj].setGradient(1.0);
+      tape.evaluate();
+      for (int ii=0; ii<state_size; ii++)
+      {
+         jac(jj,ii)=ad_state[ii].getGradient();
+      }
+      tape.clearAdjoints();
+      ad_result[jj].setGradient(0.0);
+   }
+   tape.resetTo(pos);
+#endif
+#else
+   MFEM_WARNING("not compiled with CoDiPack -- will return zero gradient");
+#endif
+}
+
+void ADVectorTDFunc::Curl(const Vector &vstate, const real_t time, Vector &curl)
+{
+   MFEM_ASSERT(vstate.Size() == 3, "state size expected to be 3");
+   MFEM_ASSERT(curl.Size() == 3, "solution size expected to be 3");
+   DenseMatrix jac(vstate.Size()); // Assume the jacobian is square
+   Jacobian(vstate, time, jac);
+
+   curl[0] = jac(2,1) - jac(1,2);
+   curl[1] = jac(0,2) - jac(2,0);
+   curl[2] = jac(1,0) - jac(0,1);
+}
+
+real_t ADVectorTDFunc::Divergence(const Vector &vstate, const real_t time)
+{
+   DenseMatrix jac(vstate.Size()); // Assume the jacobian is square
+   Jacobian(vstate, time, jac);
+   real_t div = 0.0;
+   for (int ii = 0; ii<vstate.Size(); ii++)
+   {
+      div += jac(ii,ii);
+   }
+   return div;
+}
+
+void ADVectorTDFunc::Gradient(const Vector &vstate, const real_t time,
+                              Vector &grad)
+{
+   MFEM_ASSERT( grad.Size() == vstate.Size(), "grad and state size not equal");
+
+   DenseMatrix jac(1, vstate.Size());
+   Jacobian(vstate, time, jac);
+
+   jac.GetRow(0, grad);
+}
+
+void ADVectorTDFunc::Rate(const Vector &vstate, const real_t time, Vector &rate)
+{
+#ifdef MFEM_USE_CODIPACK
+   ADVector ad_state(vstate.Size());
+   ADVector ad_result(rate.Size());
+   ADFloat  ad_t(time);
+#ifdef MFEM_USE_ADFORWARD
+   ad_t.setGradient(1.0);
+   for (int i=0; i<vstate.Size(); i++)
+   {
+      ad_state[i].setValue(vstate[i]);
+      ad_state[i].setGradient(0.0);
+   }
+
+   F(param,ad_state,ad_t,ad_result);
+   for (int jj=0; jj<rate.Size(); jj++)
+   {
+      rate[jj] = ad_result[jj].getGradient();
+   }
+#else
+   for (int i=0; i<vstate.Size(); i++)
+   {
+      ad_state[i] = vstate[i];
+   }
+
+   ADFloat::Tape& tape = ADFloat::getTape();
+   typename ADFloat::Tape::Position pos = tape.getPosition();
+
+   tape.setActive();
+   tape.registerInput(ad_t);
+   F(param,ad_state,ad_t,ad_result);
+   for (int ii=0; ii<rate.Size(); ii++) { tape.registerOutput(ad_result[ii]); }
+   tape.setPassive();
+
+   for (int jj=0; jj<rate.Size(); jj++)
+   {
+      ad_result[jj].setGradient(1.0);
+      tape.evaluate();
+
+      rate[jj] = ad_t.getGradient();
+      tape.clearAdjoints();
+      ad_result[jj].setGradient(0.0);
+   }
+   tape.resetTo(pos);
+#endif
+#else
+   MFEM_WARNING("not compiled with CoDiPack -- will return zero rate");
+   rate = 0.0;
+#endif
+}
+
+real_t ADVectorTDFunc::ScalarRate(const Vector &vstate, const real_t time)
+{
+   Vector rv(1);
+   Rate(vstate, time, rv);
+   return rv[0];
+}
+
+void ADVectorTDFunc::Solution(const Vector &vstate, const real_t time,
+                              Vector &sol)
+{
+   ADVector ad_state(vstate.Size());
+   ADVector ad_result(sol.Size());
+   ADFloat  ad_t(time);
+   for (int i=0; i<vstate.Size(); i++)
+   {
+#ifdef MFEM_USE_CODIPACK
+      ad_state[i].setValue(vstate[i]);
+      ad_state[i].setGradient(0.0);
+#else
+      ad_state[i] = vstate[i];
+#endif
+   }
+
+   F(param,ad_state,ad_t,ad_result);
+   for (int jj=0; jj<sol.Size(); jj++)
+   {
+#ifdef MFEM_USE_CODIPACK
+      sol[jj] = ad_result[jj].getValue();
+#else
+      sol[jj] = ad_result[jj];
+#endif
+   }
+}
+
+real_t ADVectorTDFunc::ScalarSolution(const Vector &vstate, const real_t time)
+{
+   Vector sol(1);
+   Solution(vstate, time, sol);
+   return sol[0];
+}
+
+}

fem/autodiff.hpp

@@ -0,0 +1,286 @@
+
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+
+#ifndef MFEM_AD
+#define MFEM_AD
+
+#include <vector>
+#include <functional>
+#include "../linalg/vector.hpp"
+#include "../linalg/densemat.hpp"
+
+#ifdef MFEM_USE_CODIPACK
+#include <codi.hpp>
+#endif
+
+using namespace std::placeholders;
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CODIPACK
+#ifdef MFEM_USE_ADFORWARD
+/// Forward AD type declaration
+typedef codi::RealForward ADFloat;
+#else
+/// Reverse AD type declaration
+typedef codi::RealReverse ADFloat;
+#endif
+#else
+/// Default to standard real
+typedef real_t ADFloat;
+#endif
+
+/// Vector type for ADFloat
+typedef std::vector<ADFloat> ADVector;
+
+/// This helper class provides evaluation of first order derivatives, namely
+/// gradient, curl and diverence, of a vector function provided in the
+/// constructor. The derivatives are evaluated with the help of automatic
+/// differentiation (AD). The derivative member functions are provided in a format
+/// compatbile with the coefficient class.
+class ADVectorFunc
+{
+private:
+   std::function<void(const Vector&, const ADVector&, ADVector&)> F;
+   Vector param;
+
+public:
+   /// F_ is user implemented function to be differentiated by
+   /// ADVectorFunc. The signature of the function is: F_(Vector&
+   /// parameters, ad::ADVectorType& state_vector, ad::ADVectorType& result).
+   /// The parameters vector should have size param_size. The state_vector
+   /// should have size state_size, and the result vector should have size
+   /// vector_size. All size parameters are teplate parameters in
+   /// ADVectorFunc.
+   ADVectorFunc(
+      std::function<void(const Vector&, const ADVector&, ADVector&)> F_,
+      Vector &par)
+   {
+      F=F_;
+      param = par;
+   }
+
+   ADVectorFunc(
+      std::function<void(const Vector&, const ADVector&, ADVector&)> F_)
+   {
+      F=F_;
+   }
+
+   /// Evaluates the Jacobian of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The function and parameters are
+   /// provided in the Constructor. The Jacobian (jac) has dimensions
+   /// [vector_size x state_size].
+   void Jacobian(const Vector &vstate,DenseMatrix &jac);
+
+   /// Evaluates the Curl of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The function and parameters are
+   /// provided in the Constructor. The Curl only works in 3D.
+   void Curl(const Vector &vstate, Vector &curl);
+
+   /// Evaluates the Divergence of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The function and parameters are
+   /// provided in the Constructor.
+   real_t Divergence(const Vector &vstate);
+
+   /// Evaluates the Gradient of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The Gradient (grad) has dimensions
+   /// [ state_size].
+   void Gradient(const Vector &vstate, Vector &grad);
+
+   /// Evaluates the Solution of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The Solution (sol) has dimensions
+   /// [ state_size]. The function and parameters are provided in the Constructor.
+   void Solution(const Vector &vstate, Vector &sol);
+
+   /// Evaluates the Solution of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate.  The function and parameters are
+   /// provided in the Constructor.
+   real_t ScalarSolution(const Vector &vstate);
+
+   /// Function type definitions
+   using ScalFunc = real_t(const Vector &);
+   using VecFunc = void(const Vector &, Vector &);
+   using MatFunc = void(const Vector &, DenseMatrix &);
+
+   /// Access function to get the Jacobian Function
+   std::function<MatFunc> GetJacobian()
+   {
+      return std::bind(&ADVectorFunc::Jacobian, *this, _1, _2);
+   }
+
+   /// Access function to get the Curl Function
+   std::function<VecFunc> GetCurl()
+   {
+      return std::bind(&ADVectorFunc::Curl, *this, _1, _2);
+   }
+
+   /// Access function to get the Divergence Function
+   std::function<ScalFunc> GetDivergence()
+   {
+      return std::bind(&ADVectorFunc::Divergence, *this, _1);
+   }
+
+   /// Access function to get the Gradient Function
+   std::function<VecFunc> GetGradient()
+   {
+      return std::bind(&ADVectorFunc::Gradient, *this, _1, _2);
+   }
+
+   /// Access function to get the Scalar Function
+   std::function<VecFunc> GetSolution()
+   {
+      return std::bind(&ADVectorFunc::Solution, *this, _1, _2);
+   }
+
+   /// Evaluates the Solution of the vector function F_ for a set of parameters
+   std::function<ScalFunc> GetScalarSolution()
+   {
+      return std::bind(&ADVectorFunc::ScalarSolution, *this, _1);
+   }
+};
+
+/// This helper class provides evaluation of first order derivatives, namely
+/// gradient, curl and diverence, of a vector function provided in the
+/// constructor. The derivatives are evaluated with the help of automatic
+/// differentiation (AD). The derivative member functions are provided in a format
+/// compatbile with the coefficient class.
+class ADVectorTDFunc
+{
+private:
+   std::function<void(const Vector&, const ADVector&, const ADFloat&, ADVector&)>
+   F;
+   Vector param;
+
+public:
+   /// F_ is user implemented function to be differentiated by
+   /// ADVectorFunc. The signature of the function is: F_(Vector&
+   /// parameters, ad::ADVectorType& state_vector, ad::ADVectorType& result).
+   /// The parameters vector should have size param_size. The state_vector
+   /// should have size state_size, and the result vector should have size
+   /// vector_size. All size parameters are teplate parameters in
+   /// ADVectorFunc.
+   ADVectorTDFunc(
+      std::function<void(const Vector&, const ADVector&, const ADFloat&, ADVector&)>
+      F_,
+      Vector &par)
+   {
+      F=F_;
+      param = par;
+   }
+
+   ADVectorTDFunc(
+      std::function<void(const Vector&, const ADVector&, const ADFloat&, ADVector&)>
+      F_)
+   {
+      F=F_;
+   }
+
+   /// Evaluates the Jacobian of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The function and parameters are
+   /// provided in the Constructor. The Jacobian (jac) has dimensions
+   /// [vector_size x state_size].
+   void Jacobian(const Vector &vstate, const real_t time, DenseMatrix &jac);
+
+   /// Evaluates the Curl of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The function and parameters are
+   /// provided in the Constructor. The Curl only works in 3D.
+   void Curl(const Vector &vstate, const real_t time, Vector &curl);
+
+   /// Evaluates the Divergence of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The function and parameters are
+   /// provided in the Constructor.
+   real_t Divergence(const Vector &vstate, const real_t time);
+
+   /// Evaluates the Gradient of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The Gradient (grad) has dimensions
+   /// [ state_size].
+   void Gradient(const Vector &vstate, const real_t time, Vector &grad);
+
+   /// Evaluates the Rate of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The Solution (sol) has dimensions
+   /// [ state_size]. The function and parameters are provided in the Constructor.
+   void Rate(const Vector &vstate, const real_t time, Vector &sol);
+
+   /// Evaluates the Rate of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate.  The function and parameters are
+   /// provided in the Constructor.
+   real_t ScalarRate(const Vector &vstate, const real_t time);
+
+   /// Evaluates the Solution of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate. The Solution (sol) has dimensions
+   /// [ state_size]. The function and parameters are provided in the Constructor.
+   void Solution(const Vector &vstate, const real_t time, Vector &sol);
+
+   /// Evaluates the Solution of the vector function F_ for a set of parameters
+   /// (vparam) and state vector vstate.  The function and parameters are
+   /// provided in the Constructor.
+   real_t ScalarSolution(const Vector &vstate, const real_t time);
+
+   /// Function type definitions
+   using ScalTDFunc = real_t(const Vector &, const real_t time);
+   using VecTDFunc = void(const Vector &, const real_t time, Vector &);
+   using MatTDFunc = void(const Vector &, const real_t time, DenseMatrix &);
+
+   /// Access function to get the Jacobian Function
+   std::function<MatTDFunc> GetJacobian()
+   {
+      return std::bind(&ADVectorTDFunc::Jacobian, *this, _1, _2, _3);
+   }
+
+   /// Access function to get the Curl Function
+   std::function<VecTDFunc> GetCurl()
+   {
+      return std::bind(&ADVectorTDFunc::Curl, *this, _1, _2, _3);
+   }
+
+   /// Access function to get the Divergence Function
+   std::function<ScalTDFunc> GetDivergence()
+   {
+      return std::bind(&ADVectorTDFunc::Divergence, *this, _1, _2);
+   }
+
+   /// Access function to get the Gradient Function
+   std::function<VecTDFunc> GetGradient()
+   {
+      return std::bind(&ADVectorTDFunc::Gradient, *this, _1, _2, _3);
+   }
+
+   /// Access function to get the Rate Function
+   std::function<VecTDFunc> GetRate()
+   {
+      return std::bind(&ADVectorTDFunc::Rate, *this, _1, _2, _3);
+   }
+
+   /// Access function to get the Scalar Rate Function
+   std::function<ScalTDFunc> GetScalarRate()
+   {
+      return std::bind(&ADVectorTDFunc::ScalarRate, *this, _1, _2);
+   }
+
+   /// Access function to get the Solution Function
+   std::function<VecTDFunc> GetSolution()
+   {
+      return std::bind(&ADVectorTDFunc::Solution, *this, _1, _2, _3);
+   }
+
+   /// Access function to get the Scalar Solution Function
+   std::function<ScalTDFunc> GetScalarSolution()
+   {
+      return std::bind(&ADVectorTDFunc::ScalarSolution, *this, _1, _2);
+   }
+
+};
+
+}
+#endif

fem/coefficient.cpp

@@ -144,6 +144,24 @@ real_t FunctionCoefficient::Eval(ElementTransformation & T,
    }
 }
 
+/*real_t ADFunctionCoefficient::Eval(ElementTransformation & T,
+                                 const IntegrationPoint & ip)
+{
+   real_t x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   if (Function)
+   {
+      return Function(transip);
+   }
+   else
+   {
+      return TDFunction(transip, GetTime());
+   }
+}*/
+
 real_t CartesianCoefficient::Eval(ElementTransformation & T,
                                   const IntegrationPoint & ip)
 {
@@ -1463,6 +1481,54 @@ real_t LpNormLoop(real_t p, VectorCoefficient &coeff, Mesh &mesh,
    return norm;
 }
 
+real_t LpNormLoop(real_t p, MatrixCoefficient &coeff, Mesh &mesh,
+                  const IntegrationRule *irs[])
+{
+   real_t norm = 0.0;
+   ElementTransformation *tr;
+   int vdim = coeff.GetVDim();
+   DenseMatrix vval(vdim);
+   real_t val;
+
+   for (int i = 0; i < mesh.GetNE(); i++)
+   {
+      tr = mesh.GetElementTransformation(i);
+      const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+         tr->SetIntPoint(&ip);
+         coeff.Eval(vval, *tr, ip);
+         if (p < infinity())
+         {
+            for (int idim(0); idim < vdim; ++idim)
+            {
+               for (int jdim(0); jdim < vdim; ++jdim)
+               {
+                  norm += ip.weight * tr->Weight() * pow(fabs( vval(idim,jdim) ), p);
+               }
+            }
+         }
+         else
+         {
+            for (int idim(0); idim < vdim; ++idim)
+            {
+               for (int jdim(0); jdim < vdim; ++jdim)
+               {
+                  val = fabs(vval(idim, jdim));
+                  if (norm < val)
+                  {
+                     norm = val;
+                  }
+               }
+            }
+         }
+      }
+   }
+
+   return norm;
+}
+
 real_t ComputeLpNorm(real_t p, Coefficient &coeff, Mesh &mesh,
                      const IntegrationRule *irs[])
 {
@@ -1505,6 +1571,27 @@ real_t ComputeLpNorm(real_t p, VectorCoefficient &coeff, Mesh &mesh,
    return norm;
 }
 
+real_t ComputeLpNorm(real_t p, MatrixCoefficient &coeff, Mesh &mesh,
+                     const IntegrationRule *irs[])
+{
+   real_t norm = LpNormLoop(p, coeff, mesh, irs);
+
+   if (p < infinity())
+   {
+      // negative quadrature weights may cause norm to be negative
+      if (norm < 0.0)
+      {
+         norm = -pow(-norm, 1.0/p);
+      }
+      else
+      {
+         norm = pow(norm, 1.0/p);
+      }
+   }
+
+   return norm;
+}
+
 #ifdef MFEM_USE_MPI
 real_t ComputeGlobalLpNorm(real_t p, Coefficient &coeff, ParMesh &pmesh,
                            const IntegrationRule *irs[])
@@ -1569,6 +1656,38 @@ real_t ComputeGlobalLpNorm(real_t p, VectorCoefficient &coeff, ParMesh &pmesh,
 
    return glob_norm;
 }
+
+real_t ComputeGlobalLpNorm(real_t p, MatrixCoefficient &coeff, ParMesh &pmesh,
+                           const IntegrationRule *irs[])
+{
+   real_t loc_norm = LpNormLoop(p, coeff, pmesh, irs);
+   real_t glob_norm = 0;
+
+   MPI_Comm comm = pmesh.GetComm();
+
+   if (p < infinity())
+   {
+      MPI_Allreduce(&loc_norm, &glob_norm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                    comm);
+
+      // negative quadrature weights may cause norm to be negative
+      if (glob_norm < 0.0)
+      {
+         glob_norm = -pow(-glob_norm, 1.0/p);
+      }
+      else
+      {
+         glob_norm = pow(glob_norm, 1.0/p);
+      }
+   }
+   else
+   {
+      MPI_Allreduce(&loc_norm, &glob_norm, 1, MPITypeMap<real_t>::mpi_type, MPI_MAX,
+                    comm);
+   }
+
+   return glob_norm;
+}
 #endif
 
 VectorQuadratureFunctionCoefficient::VectorQuadratureFunctionCoefficient(

fem/coefficient.hpp

@@ -18,6 +18,7 @@
 #include "../linalg/linalg.hpp"
 #include "intrules.hpp"
 #include "eltrans.hpp"
+#include "autodiff.hpp"
 
 namespace mfem
 {
@@ -2304,6 +2305,9 @@ inline int operator&(CoefficientStorage a, CoefficientStorage b)
 }
 
 
+
+
+
 /// @brief Class to represent a coefficient evaluated at quadrature points.
 ///
 /// In the general case, a CoefficientVector is the same as a QuadratureFunction
@@ -2406,6 +2410,11 @@ real_t ComputeLpNorm(real_t p, Coefficient &coeff, Mesh &mesh,
 real_t ComputeLpNorm(real_t p, VectorCoefficient &coeff, Mesh &mesh,
                      const IntegrationRule *irs[]);
 
+/** @brief Compute the Lp norm of a matrix function f = {f_ij}_i,j=1...N.
+    $ \| f \|_{Lp} = ( \sum_{ij} \| f_{ij} \|_{Lp}^p )^{1/p} $ */
+real_t ComputeLpNorm(real_t p, MatrixCoefficient &coeff, Mesh &mesh,
+                     const IntegrationRule *irs[]);
+
 #ifdef MFEM_USE_MPI
 /** @brief Compute the global Lp norm of a function f.
     $ \| f \|_{Lp} = ( \int_\Omega | f |^p d\Omega)^{1/p} $ */

fem/fem.hpp

@@ -20,6 +20,7 @@
 #include "doftrans.hpp"
 #include "eltrans.hpp"
 #include "coefficient.hpp"
+#include "autodiff.hpp"
 #include "complex_fem.hpp"
 #include "convergence.hpp"
 #include "lininteg.hpp"

tests/unit/CMakeLists.txt

@@ -70,6 +70,7 @@ set(UNIT_TESTS_SRCS
   fem/test_calcshape.cpp
   fem/test_calcvshape.cpp
   fem/test_coefficient.cpp
+  fem/test_autodiff.cpp
   fem/test_datacollection.cpp
   fem/test_derefine.cpp
   fem/test_dgmassinv.cpp

tests/unit/fem/test_autodiff.cpp

@@ -0,0 +1,874 @@
+// Copyright (c) 2010-2024, 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 "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+namespace autodiff
+{
+void scalarFun(const Vector& param, const ADVector& x, ADVector& f)
+{
+   int d = x.size();
+
+   real_t pi  = (real_t)(M_PI);
+   f[0] = sin(pi*x[0]);
+
+   if (d >= 2)
+   {
+      f[0] *= sin(2*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f[0] *= sin(3*pi*x[2]);
+   }
+}
+
+real_t scalarFunD(const Vector& x)
+{
+   int d = x.Size();
+
+   real_t pi  = (real_t)(M_PI);
+   real_t f = sin(pi*x[0]);
+
+   if (d >= 2)
+   {
+      f *= sin(2*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f *= sin(3*pi*x[2]);
+   }
+   return f;
+}
+
+void scalarGrad(const Vector & x, Vector & a)
+{
+   real_t pi  = (real_t)(M_PI);
+
+   real_t sx = sin(pi*x[0]);
+   real_t cx = cos(pi*x[0]);
+   real_t sy = 1.0;
+   real_t cy = 1.0;
+   real_t sz = 1.0;
+   real_t cz = 1.0;
+
+   a[0] = pi*cx;
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(2*pi*x[1]);
+      cy = cos(2*pi*x[1]);
+      a[0] = pi*cx*sy;
+      a[1] = 2*pi*sx*cy;
+   }
+   if (d >= 3)
+   {
+      sz = sin(3*pi*x[2]);
+      cz = cos(3*pi*x[2]);
+      a[0] = pi*cx*sy*sz;
+      a[1] = 2*pi*sx*cy*sz;
+      a[2] = 3*pi*sx*sy*cz;
+   }
+}
+
+void scalarFun2(const Vector& param, const ADVector& x, ADVector& f)
+{
+   int d = x.size();
+
+   real_t pi  = (real_t)(M_PI);
+   f[0] = sin(param[0]*pi*x[0]);
+
+   if (d >= 2)
+   {
+      f[0] *= sin(param[1]*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f[0] *= sin(param[2]*pi*x[2]);
+   }
+}
+
+
+TEST_CASE("Autodiff of scalar function for Coefficient",
+          "[AD Scalar]")
+{
+   Vector param(3);
+   param[0] = 1.0;
+   param[1] = 2.0;
+   param[2] = 3.0;
+   ADVectorFunc fun(scalarFun);
+   ADVectorFunc fun2(scalarFun2, param);
+
+   SECTION("1D")
+   {
+      real_t ref_sol;
+      Vector x(1), ref_grad(1), grad(1);
+
+      for (int i = 0; i < 10; i++)
+      {
+         x[0] = 0.1*i;
+
+         // Solution
+         ref_sol = scalarFunD(x);
+
+         REQUIRE(fun.ScalarSolution(x)- ref_sol == MFEM_Approx(0.0));
+         REQUIRE(fun2.ScalarSolution(x) -ref_sol == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+         // Gradient -- > Implies Jacobian
+         scalarGrad(x,ref_grad);
+
+         fun.Gradient(x,grad);
+         grad -= ref_grad;
+         REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+
+         fun2.Gradient(x,grad);
+         grad -= ref_grad;
+         REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+#endif
+      }
+   }
+
+   SECTION("2D")
+   {
+      real_t ref_sol;
+      Vector x(2), ref_grad(2), grad(2);
+
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            x[0] = 0.1*i;
+            x[1] = 0.1*j;
+
+            // Solution
+            ref_sol = scalarFunD(x);
+
+            REQUIRE(fun.ScalarSolution(x)- ref_sol == MFEM_Approx(0.0));
+            REQUIRE(fun2.ScalarSolution(x) -ref_sol == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+            // Gradient -- > Implies Jacobian
+            scalarGrad(x,ref_grad);
+
+            fun.Gradient(x,grad);
+            grad -= ref_grad;
+            REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+
+            fun2.Gradient(x,grad);
+            grad -= ref_grad;
+            REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+#endif
+         }
+      }
+   }
+
+   SECTION("3D")
+   {
+      real_t ref_sol;
+      Vector x(3), ref_grad(3), grad(3);
+
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            for (int k = 0; k < 10; k++)
+            {
+               x[0] = 0.1*i;
+               x[1] = 0.1*j;
+               x[2] = 0.1*k;
+
+               // Solution
+               ref_sol = scalarFunD(x);
+
+               REQUIRE(fun.ScalarSolution(x)- ref_sol == MFEM_Approx(0.0));
+               REQUIRE(fun2.ScalarSolution(x) -ref_sol == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+               // Gradient -- > Implies Jacobian
+               scalarGrad(x,ref_grad);
+
+               fun.Gradient(x,grad);
+               grad -= ref_grad;
+               REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+
+               fun2.Gradient(x,grad);
+               grad -= ref_grad;
+               REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+#endif
+            }
+         }
+      }
+   }
+
+   SECTION("2D Coefficient")
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(3, 3, Element::QUADRILATERAL);
+      int order_quad = 6;
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      FunctionCoefficient sol_ref_coeff(scalarFunD);
+      FunctionCoefficient sol_coeff(fun.GetScalarSolution());
+      SumCoefficient sol_dif_coeff(sol_ref_coeff, sol_coeff, -1.0);
+
+      real_t sol_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+      REQUIRE(sol_norm == MFEM_Approx(0.0));
+
+#ifdef MFEM_USE_CODIPACK
+      VectorFunctionCoefficient grad_ref_coeff(2,scalarGrad);
+      VectorFunctionCoefficient grad_coeff(2,fun.GetGradient());
+      VectorSumCoefficient grad_dif_coeff(grad_ref_coeff, grad_coeff, -1.0);
+
+      real_t grad_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+      REQUIRE(grad_norm == MFEM_Approx(0.0));
+#endif
+   }
+
+} // Test Case
+
+
+void scalarTDFun(const Vector& param, const ADVector& x, ADFloat t, ADVector& f)
+{
+   int d = x.size();
+
+   real_t pi  = (real_t)(M_PI);
+   f[0] = (2.0*t + 3.0*t*t)*sin(3*pi*x[0]);
+
+   if (d >= 2)
+   {
+      f[0] *= sin(5*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f[0] *= sin(7*pi*x[2]);
+   }
+}
+
+real_t scalarTDFunD(const Vector& x, const real_t t)
+{
+   int d = x.Size();
+
+   real_t pi  = (real_t)(M_PI);
+   real_t f= (2.0*t + 3.0*t*t)*sin(3*pi*x[0]);
+
+   if (d >= 2)
+   {
+      f *= sin(5*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f *= sin(7*pi*x[2]);
+   }
+   return f;
+}
+
+real_t  scalarTDRate(const Vector& x, const real_t t)
+{
+   int d = x.Size();
+
+   real_t pi  = (real_t)(M_PI);
+   real_t f = (2.0 + 6.0*t)*sin(3*pi*x[0]);
+
+   if (d >= 2)
+   {
+      f *= sin(5*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f *= sin(7*pi*x[2]);
+   }
+   return f;
+}
+
+void scalarTDGrad(const Vector & x, const real_t t, Vector & a)
+{
+   real_t amp = (2.0*t + 3.0*t*t);
+
+   real_t pi  = (real_t)(M_PI);
+
+   real_t sx = sin(3*pi*x[0]);
+   real_t cx = cos(3*pi*x[0]);
+   real_t sy = 1.0;
+   real_t cy = 1.0;
+   real_t sz = 1.0;
+   real_t cz = 1.0;
+
+   a[0] = amp*3*pi*cx;
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(5*pi*x[1]);
+      cy = cos(5*pi*x[1]);
+      a[0] = amp*3*pi*cx*sy;
+      a[1] = amp*5*pi*sx*cy;
+   }
+   if (d >= 3)
+   {
+      sz = sin(7*pi*x[2]);
+      cz = cos(7*pi*x[2]);
+      a[0] = amp*3*pi*cx*sy*sz;
+      a[1] = amp*5*pi*sx*cy*sz;
+      a[2] = amp*7*pi*sx*sy*cz;
+   }
+}
+
+void scalarTDFun2(const Vector& param, const ADVector& x, const ADFloat t,
+                  ADVector& f)
+{
+   ADFloat amp = (2.0*t + 3.0*t*t);
+   int d = x.size();
+
+   real_t pi  = (real_t)(M_PI);
+   f[0] = amp*sin(param[0]*pi*x[0]);
+
+   if (d >= 2)
+   {
+      f[0] *= sin(param[1]*pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      f[0] *= sin(param[2]*pi*x[2]);
+   }
+}
+
+TEST_CASE("Autodiff of TD scalar function for Coefficient",
+          "[AD TD Scalar]")
+{
+   Vector param(3);
+   param[0] = 3.0;
+   param[1] = 5.0;
+   param[2] = 7.0;
+   ADVectorTDFunc fun(scalarTDFun);
+   ADVectorTDFunc fun2(scalarTDFun2,param);
+
+   SECTION("Time dependent - 1D")
+   {
+      real_t t, ref_sol;
+      Vector x(1);
+#ifdef MFEM_USE_CODIPACK
+      real_t ref_rate;
+      Vector ref_grad(1), grad(1);
+#endif
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            t = 0.1*i;
+            x[0] = 0.1*j;
+
+            // Solution
+            ref_sol = scalarTDFunD(x,t);
+
+            REQUIRE(fun.ScalarSolution(x,t)- ref_sol  == MFEM_Approx(0.0));
+            REQUIRE(fun2.ScalarSolution(x,t) -ref_sol == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+            // Rate
+            ref_rate = scalarTDRate(x,t);
+
+            REQUIRE(fun.ScalarRate(x,t)- ref_rate  == MFEM_Approx(0.0));
+            REQUIRE(fun2.ScalarRate(x,t) -ref_rate == MFEM_Approx(0.0));
+
+            // Gradient -- > Implies Jacobian
+            scalarTDGrad(x,t,ref_grad);
+
+            fun.Gradient(x,t,grad);
+            grad -= ref_grad;
+            REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+
+            fun2.Gradient(x,t,grad);
+            grad -= ref_grad;
+            REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+#endif
+         }
+      }
+   }
+
+   SECTION("Time dependent - 2D")
+   {
+      real_t t, ref_sol;
+      Vector x(2);
+#ifdef MFEM_USE_CODIPACK
+      real_t ref_rate;
+      Vector ref_grad(2), grad(2);
+#endif
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            for (int k = 0; k < 10; k++)
+            {
+               t = 0.1*i;
+               x[0] = 0.1*j;
+               x[1] = 0.1*k;
+
+               // Solution
+               ref_sol = scalarTDFunD(x,t);
+
+               REQUIRE(fun.ScalarSolution(x,t)- ref_sol  == MFEM_Approx(0.0));
+               REQUIRE(fun2.ScalarSolution(x,t) -ref_sol == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+               // Rate
+               ref_rate = scalarTDRate(x,t);
+
+               REQUIRE(fun.ScalarRate(x,t)- ref_rate  == MFEM_Approx(0.0));
+               REQUIRE(fun2.ScalarRate(x,t) -ref_rate == MFEM_Approx(0.0));
+
+               // Gradient -- > Implies Jacobian
+               scalarTDGrad(x,t,ref_grad);
+
+               fun.Gradient(x,t,grad);
+               grad -= ref_grad;
+               REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+
+               fun2.Gradient(x,t,grad);
+               grad -= ref_grad;
+               REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+#endif
+            }
+         }
+      }
+   }
+
+   SECTION("Time dependent - 3D")
+   {
+      real_t t, ref_sol;
+      Vector x(3);
+#ifdef MFEM_USE_CODIPACK
+      real_t ref_rate;
+      Vector ref_grad(3), grad(3);
+#endif
+
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            for (int k = 0; k < 10; k++)
+            {
+               for (int l = 0; l < 10; l++)
+               {
+
+                  t = 0.1*i;
+                  x[0] = 0.1*j;
+                  x[1] = 0.1*k;
+                  x[2] = 0.1*l;
+
+                  // Solution
+                  ref_sol = scalarTDFunD(x,t);
+
+                  REQUIRE(fun.ScalarSolution(x,t)- ref_sol  == MFEM_Approx(0.0));
+                  REQUIRE(fun2.ScalarSolution(x,t) -ref_sol == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+                  // Rate
+                  ref_rate = scalarTDRate(x,t);
+
+                  REQUIRE(fun.ScalarRate(x,t)- ref_rate  == MFEM_Approx(0.0));
+                  REQUIRE(fun2.ScalarRate(x,t) -ref_rate == MFEM_Approx(0.0));
+
+                  // Gradient -- > Implies Jacobian
+                  scalarTDGrad(x,t,ref_grad);
+
+                  fun.Gradient(x,t,grad);
+                  grad -= ref_grad;
+                  REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+
+                  fun2.Gradient(x,t,grad);
+                  grad -= ref_grad;
+                  REQUIRE(grad.Norml2() == MFEM_Approx(0.0));
+#endif
+               }
+            }
+         }
+      }
+   }
+
+   SECTION("2D Coefficient")
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(3, 3, Element::QUADRILATERAL);
+      int order_quad = 6;
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      FunctionCoefficient sol_ref_coeff(scalarTDFunD);
+      FunctionCoefficient sol_coeff(fun.GetScalarSolution());
+      SumCoefficient sol_dif_coeff(sol_ref_coeff, sol_coeff, -1.0);
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         sol_ref_coeff.SetTime(t);
+         sol_coeff.SetTime(t);
+         real_t sol_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+         REQUIRE(sol_norm == MFEM_Approx(0.0));
+      }
+
+#ifdef MFEM_USE_CODIPACK
+      FunctionCoefficient rate_ref_coeff(scalarTDRate);
+      FunctionCoefficient rate_coeff(fun.GetScalarRate());
+      SumCoefficient rate_dif_coeff(rate_ref_coeff, rate_coeff, -1.0);
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         rate_ref_coeff.SetTime(t);
+         rate_coeff.SetTime(t);
+         real_t rate_norm = ComputeLpNorm(2.0, rate_dif_coeff, mesh, irs);
+         REQUIRE(rate_norm == MFEM_Approx(0.0));
+      }
+
+      VectorFunctionCoefficient grad_ref_coeff(2,scalarTDGrad);
+      VectorFunctionCoefficient grad_coeff(2,fun.GetGradient());
+      VectorSumCoefficient grad_dif_coeff(grad_ref_coeff, grad_coeff, -1.0);
+
+      real_t grad_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+      REQUIRE(grad_norm == MFEM_Approx(0.0));
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         grad_ref_coeff.SetTime(t);
+         grad_coeff.SetTime(t);
+         real_t grad_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+         REQUIRE(grad_norm == MFEM_Approx(0.0));
+      }
+#endif
+   }
+
+} // Test Case
+
+void vectorFun(const Vector& param, const ADVector& x, ADVector& f)
+{
+   real_t pi  = (real_t)(M_PI);
+   f[0] = sin(pi*x[0])*cos(pi*x[1])*3*x[2]*x[2];
+   f[1] = cos(pi*x[0])*cos(pi*x[1])*2*x[2];
+   f[2] = -pi*cos(pi*x[0])*cos(pi*x[1])*x[2]*x[2]*x[2]
+          +pi*cos(pi*x[0])*sin(pi*x[1])*x[2]*x[2];
+}
+
+void vectorFunD(const Vector& x, Vector& f)
+{
+   real_t pi  = (real_t)(M_PI);
+   f[0] = sin(pi*x[0])*cos(pi*x[1])*3*x[2]*x[2];
+   f[1] = cos(pi*x[0])*cos(pi*x[1])*2*x[2];
+   f[2] = -pi*cos(pi*x[0])*cos(pi*x[1])*x[2]*x[2]*x[2]
+          +pi*cos(pi*x[0])*sin(pi*x[1])*x[2]*x[2];
+}
+
+void vectorJac(const Vector & x, DenseMatrix & j)
+{
+   real_t pi  = (real_t)(M_PI);
+   j(0,0) = pi*cos(pi*x[0])*cos(pi*x[1])*3*x[2]*x[2];
+   j(0,1) = -pi*sin(pi*x[0])*sin(pi*x[1])*3*x[2]*x[2];
+   j(0,2) = sin(pi*x[0])*cos(pi*x[1])*6*x[2];
+
+   j(1,0) = -pi*sin(pi*x[0])*cos(pi*x[1])*2*x[2];
+   j(1,1) = -pi*cos(pi*x[0])*sin(pi*x[1])*2*x[2];
+   j(1,2) = cos(pi*x[0])*cos(pi*x[1])*2;
+
+   j(2,0) = pi*pi*sin(pi*x[0])*cos(pi*x[1])*x[2]*x[2]*x[2]
+            -pi*pi*sin(pi*x[0])*sin(pi*x[1])*x[2]*x[2];
+   j(2,1) = pi*pi*cos(pi*x[0])*sin(pi*x[1])*x[2]*x[2]*x[2]
+            +pi*pi*cos(pi*x[0])*cos(pi*x[1])*x[2]*x[2];
+   j(2,2) = -pi*cos(pi*x[0])*cos(pi*x[1])*3*x[2]*x[2]
+            +pi*cos(pi*x[0])*sin(pi*x[1])*2*x[2];
+}
+
+void vectorCurl(const Vector & x, Vector& curl)
+{
+   DenseMatrix jac(x.Size());
+   vectorJac(x,jac);
+
+   // Curl
+   curl[0] = jac(2,1) - jac(1,2);
+   curl[1] = jac(0,2) - jac(2,0);
+   curl[2] = jac(1,0) - jac(0,1);
+}
+
+TEST_CASE("Autodiff of vector function for Coefficient",
+          "[AD Vector]")
+{
+
+   ADVectorFunc fun(vectorFun);
+
+   SECTION("3D")
+   {
+      Vector x(3), ref_sol(3), sol(3);
+#ifdef MFEM_USE_CODIPACK
+      real_t ref_div, div;
+      Vector  ref_curl(3), curl(3);
+      DenseMatrix ref_jac(3,3), jac(3,3);
+#endif
+
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            for (int k = 0; k < 10; k++)
+            {
+               x[0] = 0.1*i;
+               x[1] = 0.1*j;
+               x[2] = 0.1*k;
+
+               // Solution
+               vectorFunD(x,ref_sol);
+
+               fun.Solution(x,sol);
+               sol -= ref_sol;
+               REQUIRE(sol.Norml2() == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+               // Jacobian
+               vectorJac(x,ref_jac);
+
+               fun.Jacobian(x,jac);
+               jac -= ref_jac;
+               REQUIRE(jac.FNorm() == MFEM_Approx(0.0));
+
+               // Curl
+               ref_curl[0] = ref_jac(2,1) - ref_jac(1,2);
+               ref_curl[1] = ref_jac(0,2) - ref_jac(2,0);
+               ref_curl[2] = ref_jac(1,0) - ref_jac(0,1);
+               fun.Curl(x,curl);
+               curl -= ref_curl;
+               // REQUIRE(curl.Norml2() == MFEM_Approx(0.0));
+
+               // Divergence
+               ref_div = 0.0;
+
+               div = fun.Divergence(x);
+               REQUIRE(div-ref_div == MFEM_Approx(0.0));
+#endif
+            }
+         }
+      }
+   }
+
+   SECTION("Coefficient")
+   {
+      Mesh mesh = Mesh::MakeCartesian3D(2, 2, 2,Element::HEXAHEDRON);
+      int order_quad = 6;
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      VectorFunctionCoefficient sol_ref_coeff(3,vectorFunD);
+      VectorFunctionCoefficient sol_coeff(3,fun.GetSolution());
+      VectorSumCoefficient sol_dif_coeff(sol_ref_coeff, sol_coeff, -1.0);
+
+      real_t sol_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+      REQUIRE(sol_norm == MFEM_Approx(0.0));
+
+#ifdef MFEM_USE_CODIPACK
+      FunctionCoefficient div_coeff(fun.GetDivergence());
+
+      real_t rate_norm = ComputeLpNorm(2.0, div_coeff, mesh, irs);
+      REQUIRE(rate_norm == MFEM_Approx(0.0));
+
+      VectorFunctionCoefficient curl_ref_coeff(3,vectorCurl);
+      VectorFunctionCoefficient curl_coeff(3,fun.GetCurl());
+      VectorSumCoefficient curl_dif_coeff(curl_ref_coeff, curl_coeff, -1.0);
+
+      real_t curl_norm = ComputeLpNorm(2.0, curl_dif_coeff, mesh, irs);
+      REQUIRE(curl_norm == MFEM_Approx(0.0));
+
+      MatrixFunctionCoefficient jac_ref_coeff(3,vectorJac);
+      MatrixFunctionCoefficient jac_coeff(3,fun.GetJacobian());
+      MatrixSumCoefficient jac_dif_coeff(jac_ref_coeff, jac_coeff, -1.0);
+
+      real_t jac_norm = ComputeLpNorm(2.0, jac_dif_coeff, mesh, irs);
+      REQUIRE(jac_norm == MFEM_Approx(0.0));
+#endif
+
+   }
+} // Test Case
+
+void vectorTDFun(const Vector& param, const ADVector& x, const ADFloat t,
+                 ADVector& f)
+{
+   f[0] = x[1] * exp(2*t);
+   f[1] = x[0] * exp(2*t);
+   f[2] = x[2] * t;
+}
+
+void vectorTDFunD(const Vector& x, const real_t t, Vector& f)
+{
+   f[0] = x[1] * exp(2*t);
+   f[1] = x[0] * exp(2*t);
+   f[2] = x[2] * t;
+}
+
+void vectorTDRate(const Vector& x, const real_t t, Vector& r)
+{
+   r[0] = x[1] * 2*exp(2*t);
+   r[1] = x[0] * 2*exp(2*t);
+   r[2] = x[2];
+}
+
+void vectorTDJac(const Vector & x, const real_t t, DenseMatrix & j)
+{
+   real_t amp = exp(2*t);
+   j(0,0) = 0.0;
+   j(0,1) = amp;
+   j(0,2) = 0.0;
+
+   j(1,0) = amp;
+   j(1,1) = 0.0;
+   j(1,2) = 0.0;
+
+   j(2,0) = 0.0;
+   j(2,1) = 0.0;
+   j(2,2) = t;
+}
+
+void vectorTDCurl(const Vector & x, const real_t t, Vector& curl)
+{
+   DenseMatrix jac(x.Size());
+   vectorTDJac(x,t,jac);
+
+   // Curl
+   curl[0] = jac(2,1) - jac(1,2);
+   curl[1] = jac(0,2) - jac(2,0);
+   curl[2] = jac(1,0) - jac(0,1);
+}
+
+TEST_CASE("Autodiff of TD vector function for Coefficient",
+          "[AD TD Vector]")
+{
+   ADVectorTDFunc fun(vectorTDFun);
+
+   SECTION("3D")
+   {
+      real_t t;
+      Vector x(3), ref_sol(3), sol(3);
+#ifdef MFEM_USE_CODIPACK
+      real_t ref_div, div;
+      Vector ref_rate(3), rate(3), ref_curl(3), curl(3);
+      DenseMatrix ref_jac(3,3), jac(3,3);
+#endif
+      for (int i = 0; i < 10; i++)
+      {
+         for (int j = 0; j < 10; j++)
+         {
+            for (int k = 0; k < 10; k++)
+            {
+               for (int l = 0; l < 10; l++)
+               {
+
+                  t = 0.1*i;
+                  x[0] = 0.1*j;
+                  x[1] = 0.1*k;
+                  x[2] = 0.1*l;
+
+                  // Solution
+                  vectorTDFunD(x,t,ref_sol);
+
+                  fun.Solution(x,t,sol);
+                  sol -= ref_sol;
+                  REQUIRE(sol.Norml2() == MFEM_Approx(0.0));
+#ifdef MFEM_USE_CODIPACK
+                  // Rate
+                  vectorTDRate(x,t,ref_rate);
+
+                  fun.Rate(x,t,rate);
+                  rate -= ref_rate;
+                  REQUIRE(rate.Norml2() == MFEM_Approx(0.0));
+
+                  // Jacobian
+                  vectorTDJac(x,t,ref_jac);
+
+                  fun.Jacobian(x,t,jac);
+                  jac -= ref_jac;
+                  REQUIRE(jac.FNorm() == MFEM_Approx(0.0));
+
+                  // Curl
+                  ref_curl = 0.0;
+                  fun.Curl(x,t,curl);
+                  curl -= ref_curl;
+                  REQUIRE(curl.Norml2() == MFEM_Approx(0.0));
+
+                  // Divergence
+                  ref_div = ref_jac(0,0) + ref_jac(1,1) + ref_jac(2,2);
+
+                  div = fun.Divergence(x,t);
+                  REQUIRE(div-ref_div == MFEM_Approx(0.0));
+#endif
+               }
+            }
+         }
+      }
+   }
+
+   SECTION("Coefficient")
+   {
+      Mesh mesh = Mesh::MakeCartesian3D(2, 2, 2,Element::HEXAHEDRON);
+      int order_quad = 6;
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      VectorFunctionCoefficient sol_ref_coeff(3,vectorTDFunD);
+      VectorFunctionCoefficient sol_coeff(3,fun.GetSolution());
+      VectorSumCoefficient sol_dif_coeff(sol_ref_coeff, sol_coeff, -1.0);
+
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         sol_ref_coeff.SetTime(t);
+         sol_coeff.SetTime(t);
+         real_t sol_norm = ComputeLpNorm(2.0, sol_dif_coeff, mesh, irs);
+         REQUIRE(sol_norm == MFEM_Approx(0.0));
+      }
+
+#ifdef MFEM_USE_CODIPACK
+      FunctionCoefficient div_coeff(fun.GetDivergence());
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         sol_ref_coeff.SetTime(t);
+         sol_coeff.SetTime(t);
+         real_t rate_norm = ComputeLpNorm(2.0, div_coeff, mesh, irs);
+         REQUIRE(rate_norm == MFEM_Approx(0.0));
+      }
+
+      VectorFunctionCoefficient curl_ref_coeff(3,vectorTDCurl);
+      VectorFunctionCoefficient curl_coeff(3,fun.GetCurl());
+      VectorSumCoefficient curl_dif_coeff(curl_ref_coeff, curl_coeff, -1.0);
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         curl_ref_coeff.SetTime(t);
+         curl_coeff.SetTime(t);
+         real_t curl_norm = ComputeLpNorm(2.0, curl_dif_coeff, mesh, irs);
+         REQUIRE(curl_norm == MFEM_Approx(0.0));
+      }
+
+      MatrixFunctionCoefficient jac_ref_coeff(3,vectorTDJac);
+      MatrixFunctionCoefficient jac_coeff(3,fun.GetJacobian());
+      MatrixSumCoefficient jac_dif_coeff(jac_ref_coeff, jac_coeff, -1.0);
+      for (int i = 0; i < 10; i++)
+      {
+         real_t t = 0.1*i;
+         jac_ref_coeff.SetTime(t);
+         jac_coeff.SetTime(t);
+         real_t jac_norm = ComputeLpNorm(2.0, jac_dif_coeff, mesh, irs);
+         REQUIRE(jac_norm == MFEM_Approx(0.0));
+      }
+#endif
+   }
+} // Test Case
+
+} // namespace autodiff
+