diff --git a/framework/include/base/Assembly.h b/framework/include/base/Assembly.h
index ceed67d1a852..2082192cf765 100644
--- a/framework/include/base/Assembly.h
+++ b/framework/include/base/Assembly.h
@@ -1328,6 +1328,10 @@ class Assembly
   {
     return _vector_curl_phi;
   }
+  const VectorVariablePhiDivergence & divPhi(const MooseVariableField<RealVectorValue> &) const
+  {
+    return _vector_div_phi;
+  }
 
   const VectorVariablePhiValue & phiFace(const MooseVariableField<RealVectorValue> &) const
   {
@@ -1345,6 +1349,10 @@ class Assembly
   {
     return _vector_curl_phi_face;
   }
+  const VectorVariablePhiDivergence & divPhiFace(const MooseVariableField<RealVectorValue> &) const
+  {
+    return _vector_div_phi_face;
+  }
 
   const VectorVariablePhiValue & phiNeighbor(const MooseVariableField<RealVectorValue> &) const
   {
@@ -1364,6 +1372,10 @@ class Assembly
   {
     return _vector_curl_phi_neighbor;
   }
+  const VectorVariablePhiDivergence & divPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
+  {
+    return _vector_div_phi_neighbor;
+  }
 
   const VectorVariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
   {
@@ -1384,6 +1396,11 @@ class Assembly
   {
     return _vector_curl_phi_face_neighbor;
   }
+  const VectorVariablePhiDivergence &
+  divPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
+  {
+    return _vector_div_phi_face_neighbor;
+  }
 
   // Writeable references
   VariablePhiValue & phi(const MooseVariableField<Real> &) { return _phi; }
@@ -1431,6 +1448,10 @@ class Assembly
   {
     return _vector_curl_phi;
   }
+  VectorVariablePhiDivergence & divPhi(const MooseVariableField<RealVectorValue> &)
+  {
+    return _vector_div_phi;
+  }
 
   VectorVariablePhiValue & phiFace(const MooseVariableField<RealVectorValue> &)
   {
@@ -1448,6 +1469,10 @@ class Assembly
   {
     return _vector_curl_phi_face;
   }
+  VectorVariablePhiDivergence & divPhiFace(const MooseVariableField<RealVectorValue> &)
+  {
+    return _vector_div_phi_face;
+  }
 
   VectorVariablePhiValue & phiNeighbor(const MooseVariableField<RealVectorValue> &)
   {
@@ -1465,6 +1490,10 @@ class Assembly
   {
     return _vector_curl_phi_neighbor;
   }
+  VectorVariablePhiDivergence & divPhiNeighbor(const MooseVariableField<RealVectorValue> &)
+  {
+    return _vector_div_phi_neighbor;
+  }
   VectorVariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
   {
     return _vector_phi_face_neighbor;
@@ -1481,6 +1510,10 @@ class Assembly
   {
     return _vector_curl_phi_face_neighbor;
   }
+  VectorVariablePhiDivergence & divPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
+  {
+    return _vector_div_phi_face_neighbor;
+  }
 
   // Writeable references with array variable
   VariablePhiValue & phi(const MooseVariableField<RealEigenVector> &) { return _phi; }
@@ -1670,6 +1703,38 @@ class Assembly
     return _fe_shape_data_face_neighbor[type]->_curl_phi;
   }
 
+  template <typename OutputType>
+  const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhi(FEType type) const
+  {
+    _need_div[type] = true;
+    buildFE(type);
+    return _fe_shape_data[type]->_div_phi;
+  }
+
+  template <typename OutputType>
+  const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhiFace(FEType type) const
+  {
+    _need_div[type] = true;
+    buildFaceFE(type);
+    return _fe_shape_data_face[type]->_div_phi;
+  }
+
+  template <typename OutputType>
+  const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhiNeighbor(FEType type) const
+  {
+    _need_div[type] = true;
+    buildNeighborFE(type);
+    return _fe_shape_data_neighbor[type]->_div_phi;
+  }
+
+  template <typename OutputType>
+  const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhiFaceNeighbor(FEType type) const
+  {
+    _need_div[type] = true;
+    buildFaceNeighborFE(type);
+    return _fe_shape_data_face_neighbor[type]->_div_phi;
+  }
+
   /// On-demand computation of volume element accounting for RZ/RSpherical
   Real elementVolume(const Elem * elem) const;
 
@@ -2597,21 +2662,25 @@ class Assembly
   VectorVariablePhiGradient _vector_grad_phi;
   VectorVariablePhiSecond _vector_second_phi;
   VectorVariablePhiCurl _vector_curl_phi;
+  VectorVariablePhiDivergence _vector_div_phi;
 
   VectorVariablePhiValue _vector_phi_face;
   VectorVariablePhiGradient _vector_grad_phi_face;
   VectorVariablePhiSecond _vector_second_phi_face;
   VectorVariablePhiCurl _vector_curl_phi_face;
+  VectorVariablePhiDivergence _vector_div_phi_face;
 
   VectorVariablePhiValue _vector_phi_neighbor;
   VectorVariablePhiGradient _vector_grad_phi_neighbor;
   VectorVariablePhiSecond _vector_second_phi_neighbor;
   VectorVariablePhiCurl _vector_curl_phi_neighbor;
+  VectorVariablePhiDivergence _vector_div_phi_neighbor;
 
   VectorVariablePhiValue _vector_phi_face_neighbor;
   VectorVariablePhiGradient _vector_grad_phi_face_neighbor;
   VectorVariablePhiSecond _vector_second_phi_face_neighbor;
   VectorVariablePhiCurl _vector_curl_phi_face_neighbor;
+  VectorVariablePhiDivergence _vector_div_phi_face_neighbor;
 
   class FEShapeData
   {
@@ -2620,6 +2689,7 @@ class Assembly
     VariablePhiGradient _grad_phi;
     VariablePhiSecond _second_phi;
     VariablePhiCurl _curl_phi;
+    VariablePhiDivergence _div_phi;
   };
 
   class VectorFEShapeData
@@ -2629,6 +2699,7 @@ class Assembly
     VectorVariablePhiGradient _grad_phi;
     VectorVariablePhiSecond _second_phi;
     VectorVariablePhiCurl _curl_phi;
+    VectorVariablePhiDivergence _div_phi;
   };
 
   /// Shape function values, gradients, second derivatives for each FE type
@@ -2735,6 +2806,7 @@ class Assembly
   mutable std::map<FEType, bool> _need_second_derivative;
   mutable std::map<FEType, bool> _need_second_derivative_neighbor;
   mutable std::map<FEType, bool> _need_curl;
+  mutable std::map<FEType, bool> _need_div;
 
   /// The map from global index to variable scaling factor
   const NumericVector<Real> * _scaling_vector = nullptr;
@@ -2893,6 +2965,22 @@ template <>
 const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
 Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
 
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhi<VectorValue<Real>>(FEType type) const;
+
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const;
+
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const;
+
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
+
 template <>
 inline const ADTemplateVariablePhiGradient<RealVectorValue> &
 Assembly::adGradPhi<RealVectorValue>(const MooseVariableFE<RealVectorValue> & v) const
diff --git a/framework/include/utils/MooseTypes.h b/framework/include/utils/MooseTypes.h
index 97e4720bc368..25ccd30cd4e8 100644
--- a/framework/include/utils/MooseTypes.h
+++ b/framework/include/utils/MooseTypes.h
@@ -301,43 +301,52 @@ typedef typename OutputTools<Real>::VariableValue VariableValue;
 typedef typename OutputTools<Real>::VariableGradient VariableGradient;
 typedef typename OutputTools<Real>::VariableSecond VariableSecond;
 typedef typename OutputTools<Real>::VariableCurl VariableCurl;
+typedef typename OutputTools<Real>::VariableDivergence VariableDivergence;
 typedef typename OutputTools<Real>::VariablePhiValue VariablePhiValue;
 typedef typename OutputTools<Real>::VariablePhiGradient VariablePhiGradient;
 typedef typename OutputTools<Real>::VariablePhiSecond VariablePhiSecond;
 typedef typename OutputTools<Real>::VariablePhiCurl VariablePhiCurl;
+typedef typename OutputTools<Real>::VariablePhiDivergence VariablePhiDivergence;
 typedef typename OutputTools<Real>::VariableTestValue VariableTestValue;
 typedef typename OutputTools<Real>::VariableTestGradient VariableTestGradient;
 typedef typename OutputTools<Real>::VariableTestSecond VariableTestSecond;
 typedef typename OutputTools<Real>::VariableTestCurl VariableTestCurl;
+typedef typename OutputTools<Real>::VariableTestDivergence VariableTestDivergence;
 
 // types for vector variable
 typedef typename OutputTools<RealVectorValue>::VariableValue VectorVariableValue;
 typedef typename OutputTools<RealVectorValue>::VariableGradient VectorVariableGradient;
 typedef typename OutputTools<RealVectorValue>::VariableSecond VectorVariableSecond;
 typedef typename OutputTools<RealVectorValue>::VariableCurl VectorVariableCurl;
+typedef typename OutputTools<RealVectorValue>::VariableDivergence VectorVariableDivergence;
 typedef typename OutputTools<RealVectorValue>::VariablePhiValue VectorVariablePhiValue;
 typedef typename OutputTools<RealVectorValue>::VariablePhiGradient VectorVariablePhiGradient;
 typedef typename OutputTools<RealVectorValue>::VariablePhiSecond VectorVariablePhiSecond;
 typedef typename OutputTools<RealVectorValue>::VariablePhiCurl VectorVariablePhiCurl;
+typedef typename OutputTools<RealVectorValue>::VariablePhiDivergence VectorVariablePhiDivergence;
 typedef typename OutputTools<RealVectorValue>::VariableTestValue VectorVariableTestValue;
 typedef typename OutputTools<RealVectorValue>::VariableTestGradient VectorVariableTestGradient;
 typedef typename OutputTools<RealVectorValue>::VariableTestSecond VectorVariableTestSecond;
 typedef typename OutputTools<RealVectorValue>::VariableTestCurl VectorVariableTestCurl;
+typedef typename OutputTools<RealVectorValue>::VariableTestDivergence VectorVariableTestDivergence;
 
 // types for array variable
 typedef typename OutputTools<RealEigenVector>::VariableValue ArrayVariableValue;
 typedef typename OutputTools<RealEigenVector>::VariableGradient ArrayVariableGradient;
 typedef typename OutputTools<RealEigenVector>::VariableSecond ArrayVariableSecond;
 typedef typename OutputTools<RealEigenVector>::VariableCurl ArrayVariableCurl;
+typedef typename OutputTools<RealEigenVector>::VariableDivergence ArrayVariableDivergence;
 typedef typename OutputTools<RealEigenVector>::VariablePhiValue ArrayVariablePhiValue;
 typedef typename OutputTools<RealEigenVector>::VariablePhiGradient ArrayVariablePhiGradient;
 typedef std::vector<std::vector<Eigen::Map<RealDIMValue>>> MappedArrayVariablePhiGradient;
 typedef typename OutputTools<RealEigenVector>::VariablePhiSecond ArrayVariablePhiSecond;
 typedef typename OutputTools<RealEigenVector>::VariablePhiCurl ArrayVariablePhiCurl;
+typedef typename OutputTools<RealEigenVector>::VariablePhiDivergence ArrayVariablePhiDivergence;
 typedef typename OutputTools<RealEigenVector>::VariableTestValue ArrayVariableTestValue;
 typedef typename OutputTools<RealEigenVector>::VariableTestGradient ArrayVariableTestGradient;
 typedef typename OutputTools<RealEigenVector>::VariableTestSecond ArrayVariableTestSecond;
 typedef typename OutputTools<RealEigenVector>::VariableTestCurl ArrayVariableTestCurl;
+typedef typename OutputTools<RealEigenVector>::VariableTestDivergence ArrayVariableTestDivergence;
 
 /**
  * AD typedefs
diff --git a/framework/include/variables/MooseVariableData.h b/framework/include/variables/MooseVariableData.h
index 24439d5baaba..f419c0bce684 100644
--- a/framework/include/variables/MooseVariableData.h
+++ b/framework/include/variables/MooseVariableData.h
@@ -190,6 +190,16 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
    */
   const FieldVariablePhiCurl & curlPhiFace() const;
 
+  /**
+   * divergence_phi getter
+   */
+  const FieldVariablePhiDivergence & divPhi() const;
+
+  /**
+   * divergence_phi_face getter
+   */
+  const FieldVariablePhiDivergence & divPhiFace() const;
+
   /**
    * ad_grad_phi getter
    */
@@ -223,6 +233,11 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
    */
   bool computingCurl() const { return _need_curl || _need_curl_old; }
 
+  /**
+   * Whether or not this variable is computing the divergence
+   */
+  bool computingDiv() const { return _need_div || _need_div_old; }
+
   //////////////////////////////// Nodal stuff ///////////////////////////////////////////
 
   bool isNodal() const override { return _is_nodal; }
@@ -271,6 +286,12 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
    */
   const FieldVariableCurl & curlSln(Moose::SolutionState state) const;
 
+  /**
+   * Local solution divergence getter
+   * @param state The state of the simulation: current, old, older
+   */
+  const FieldVariableDivergence & divSln(Moose::SolutionState state) const;
+
   const ADTemplateVariableValue<OutputType> & adSln() const
   {
     _need_ad = _need_ad_u = true;
@@ -476,6 +497,11 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
   mutable bool _need_curl_old;
   mutable bool _need_curl_older;
 
+  /// divergence flags
+  mutable bool _need_div;
+  mutable bool _need_div_old;
+  mutable bool _need_div_older;
+
   /// AD flags
   mutable bool _need_ad;
   mutable bool _need_ad_u;
@@ -500,6 +526,11 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
   FieldVariableCurl _curl_u_old;
   FieldVariableCurl _curl_u_older;
 
+  /// divergence_u
+  FieldVariableDivergence _div_u;
+  FieldVariableDivergence _div_u_old;
+  FieldVariableDivergence _div_u_older;
+
   /// AD u
   ADTemplateVariableValue<OutputType> _ad_u;
   ADTemplateVariableGradient<OutputType> _ad_grad_u;
@@ -542,6 +573,7 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
   const FieldVariablePhiGradient * _grad_phi;
   mutable const FieldVariablePhiSecond * _second_phi;
   mutable const FieldVariablePhiCurl * _curl_phi;
+  mutable const FieldVariablePhiDivergence * _div_phi;
 
   // Mapped array phi
   MappedArrayVariablePhiGradient _mapped_grad_phi;
@@ -554,6 +586,7 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
   const FieldVariablePhiGradient * _grad_phi_face;
   mutable const FieldVariablePhiSecond * _second_phi_face;
   mutable const FieldVariablePhiCurl * _curl_phi_face;
+  mutable const FieldVariablePhiDivergence * _div_phi_face;
 
   const ADTemplateVariablePhiGradient<OutputShape> * _ad_grad_phi;
   const ADTemplateVariablePhiGradient<OutputShape> * _ad_grad_phi_face;
@@ -563,6 +596,7 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
   const FieldVariablePhiGradient * _current_grad_phi;
   const FieldVariablePhiSecond * _current_second_phi;
   const FieldVariablePhiCurl * _current_curl_phi;
+  const FieldVariablePhiDivergence * _current_div_phi;
   const ADTemplateVariablePhiGradient<OutputShape> * _current_ad_grad_phi;
 
   // dual mortar
@@ -596,6 +630,14 @@ class MooseVariableData : public MooseVariableDataBase<OutputType>
                                                                            FEType)>
       _curl_phi_face_assembly_method;
 
+  std::function<const typename OutputTools<OutputShape>::VariablePhiDivergence &(const Assembly &,
+                                                                           FEType)>
+      _div_phi_assembly_method;
+
+  std::function<const typename OutputTools<OutputShape>::VariablePhiDivergence &(const Assembly &,
+                                                                           FEType)>
+      _div_phi_face_assembly_method;
+
   std::function<const ADTemplateVariablePhiGradient<OutputShape> &(const Assembly &, FEType)>
       _ad_grad_phi_assembly_method;
   std::function<const ADTemplateVariablePhiGradient<OutputShape> &(const Assembly &, FEType)>
diff --git a/framework/include/variables/MooseVariableFE.h b/framework/include/variables/MooseVariableFE.h
index f8d232d3b150..f78f969d6fb6 100644
--- a/framework/include/variables/MooseVariableFE.h
+++ b/framework/include/variables/MooseVariableFE.h
@@ -140,6 +140,11 @@ class MooseVariableFE : public MooseVariableField<OutputType>
    */
   bool computingCurl() const override final;
 
+  /**
+   * Whether or not this variable is computing the divergence
+   */
+  bool computingDiv() const override final;
+
   bool isNodal() const override { return _element_data->isNodal(); }
   bool hasDoFsOnNodes() const override { return _element_data->hasDoFsOnNodes(); }
   FEContinuity getContinuity() const override { return _element_data->getContinuity(); };
@@ -197,6 +202,7 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   }
   const FieldVariablePhiSecond & secondPhi() const override final;
   const FieldVariablePhiCurl & curlPhi() const override final;
+  const FieldVariablePhiDivergence & divPhi() const override final;
 
   const FieldVariablePhiValue & phiFace() const override final { return _element_data->phiFace(); }
   const FieldVariablePhiGradient & gradPhiFace() const override final
@@ -209,6 +215,7 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   }
   const FieldVariablePhiSecond & secondPhiFace() const override final;
   const FieldVariablePhiCurl & curlPhiFace() const;
+  const FieldVariablePhiDivergence & divPhiFace() const;
 
   const FieldVariablePhiValue & phiNeighbor() const override final { return _neighbor_data->phi(); }
   const FieldVariablePhiGradient & gradPhiNeighbor() const override final
@@ -221,6 +228,7 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   }
   const FieldVariablePhiSecond & secondPhiNeighbor() const override final;
   const FieldVariablePhiCurl & curlPhiNeighbor() const;
+  const FieldVariablePhiDivergence & divPhiNeighbor() const;
 
   const FieldVariablePhiValue & phiFaceNeighbor() const override final
   {
@@ -236,6 +244,7 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   }
   const FieldVariablePhiSecond & secondPhiFaceNeighbor() const override final;
   const FieldVariablePhiCurl & curlPhiFaceNeighbor() const;
+  const FieldVariablePhiDivergence & divPhiFaceNeighbor() const;
 
   virtual const FieldVariablePhiValue & phiLower() const override { return _lower_data->phi(); }
   const FieldVariablePhiGradient & gradPhiLower() const { return _lower_data->gradPhi(); }
@@ -320,6 +329,11 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   const FieldVariableCurl & curlSlnOld() const { return _element_data->curlSln(Moose::Old); }
   const FieldVariableCurl & curlSlnOlder() const { return _element_data->curlSln(Moose::Older); }
 
+  /// element divergence
+  const FieldVariableDivergence & divSln() const { return _element_data->divSln(Moose::Current); }
+  const FieldVariableDivergence & divSlnOld() const { return _element_data->divSln(Moose::Old); }
+  const FieldVariableDivergence & divSlnOlder() const { return _element_data->divSln(Moose::Older); }
+
   /// AD
   const ADTemplateVariableValue<OutputType> & adSln() const override
   {
@@ -453,6 +467,20 @@ class MooseVariableFE : public MooseVariableField<OutputType>
     return _neighbor_data->curlSln(Moose::Older);
   }
 
+  /// neighbor solution divergence
+  const FieldVariableDivergence & divSlnNeighbor() const
+  {
+    return _neighbor_data->divSln(Moose::Current);
+  }
+  const FieldVariableDivergence & divSlnOldNeighbor() const
+  {
+    return _neighbor_data->divSln(Moose::Old);
+  }
+  const FieldVariableDivergence & divSlnOlderNeighbor() const
+  {
+    return _neighbor_data->divSln(Moose::Older);
+  }
+
   /// neighbor dots
   const FieldVariableValue & uDotNeighbor() const { return _neighbor_data->uDot(); }
   const FieldVariableValue & uDotDotNeighbor() const { return _neighbor_data->uDotDot(); }
diff --git a/framework/include/variables/MooseVariableFV.h b/framework/include/variables/MooseVariableFV.h
index 259fbc4bec51..5575f9741e71 100644
--- a/framework/include/variables/MooseVariableFV.h
+++ b/framework/include/variables/MooseVariableFV.h
@@ -71,6 +71,7 @@ class MooseVariableFV : public MooseVariableField<OutputType>
   using DoFValue = typename MooseVariableField<OutputType>::DoFValue;
 
   using FieldVariablePhiValue = typename MooseVariableField<OutputType>::FieldVariablePhiValue;
+  using FieldVariablePhiDivergence = typename MooseVariableField<OutputType>::FieldVariablePhiDivergence;
   using FieldVariablePhiGradient =
       typename MooseVariableField<OutputType>::FieldVariablePhiGradient;
   using FieldVariablePhiSecond = typename MooseVariableField<OutputType>::FieldVariablePhiSecond;
@@ -565,6 +566,7 @@ class MooseVariableFV : public MooseVariableField<OutputType>
 
   bool computingSecond() const override final { return false; }
   bool computingCurl() const override final { return false; }
+  bool computingDiv() const override final { return false; }
   bool usesSecondPhiNeighbor() const override final { return false; }
 
   const FieldVariablePhiValue & phi() const override final { return _phi; }
@@ -577,6 +579,10 @@ class MooseVariableFV : public MooseVariableField<OutputType>
   {
     mooseError("We don't currently implement curl for FV");
   }
+  const FieldVariablePhiDivergence & divPhi() const override final
+  {
+    mooseError("We don't currently implement divergence for FV");
+  }
 
   const FieldVariablePhiValue & phiFace() const override final { return _phi_face; }
   const FieldVariablePhiGradient & gradPhiFace() const override final { return _grad_phi_face; }
diff --git a/framework/include/variables/MooseVariableField.h b/framework/include/variables/MooseVariableField.h
index 5480fc7c255e..b378ea84a4f9 100644
--- a/framework/include/variables/MooseVariableField.h
+++ b/framework/include/variables/MooseVariableField.h
@@ -216,6 +216,11 @@ class MooseVariableField : public MooseVariableFieldBase,
    */
   virtual bool computingCurl() const = 0;
 
+    /**
+   * Whether or not this variable is computing any divergence quantities
+   */
+  virtual bool computingDiv() const = 0;
+
   /**
    * Return the variable's elemental shape functions
    */
@@ -236,6 +241,11 @@ class MooseVariableField : public MooseVariableFieldBase,
    */
   virtual const FieldVariablePhiValue & curlPhi() const = 0;
 
+  /**
+   * Divergence of the shape functions
+   */
+  virtual const FieldVariablePhiDivergence & divPhi() const = 0;
+
   /**
    * Return the variable's shape functions on an element face
    */
diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C
index 99a85ba324a7..8c991286bc6f 100644
--- a/framework/src/base/Assembly.C
+++ b/framework/src/base/Assembly.C
@@ -796,6 +796,9 @@ Assembly::reinitFE(const Elem * elem)
     if (_need_curl.find(fe_type) != _need_curl.end())
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_curl_phi()));
+    if (_need_div.find(fe_type) != _need_div.end())
+      fesd._div_phi.shallowCopy(
+          const_cast<std::vector<std::vector<Real>> &>(fe.get_div_phi()));
   }
   if (!_unique_fe_helper.empty())
   {
@@ -1296,6 +1299,9 @@ Assembly::reinitFEFace(const Elem * elem, unsigned int side)
     if (_need_curl.find(fe_type) != _need_curl.end())
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
+    if (_need_div.find(fe_type) != _need_div.end())
+      fesd._div_phi.shallowCopy(
+          const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
   }
   if (!_unique_fe_face_helper.empty())
   {
@@ -1596,6 +1602,9 @@ Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> &
     if (_need_curl.find(fe_type) != _need_curl.end())
       fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
           fe_face_neighbor.get_curl_phi()));
+    if (_need_div.find(fe_type) != _need_div.end())
+      fesd._div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(
+          fe_face_neighbor.get_div_phi()));
   }
   if (!_unique_fe_face_neighbor_helper.empty())
   {
@@ -1649,6 +1658,9 @@ Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & ref
     if (_need_curl.find(fe_type) != _need_curl.end())
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_curl_phi()));
+    if (_need_div.find(fe_type) != _need_div.end())
+      fesd._div_phi.shallowCopy(
+          const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_div_phi()));
   }
   if (!_unique_fe_neighbor_helper.empty())
   {
@@ -2055,6 +2067,9 @@ Assembly::reinitElemFaceRef(const Elem * elem,
     if (_need_curl.find(fe_type) != _need_curl.end())
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
+    if (_need_div.find(fe_type) != _need_div.end())
+      fesd._div_phi.shallowCopy(
+          const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
   }
   if (!_unique_fe_face_helper.empty())
   {
@@ -2219,6 +2234,9 @@ Assembly::reinitNeighborFaceRef(const Elem * neighbor,
     if (_need_curl.find(fe_type) != _need_curl.end())
       fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
           fe_face_neighbor.get_curl_phi()));
+    if (_need_div.find(fe_type) != _need_div.end())
+      fesd._div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(
+          fe_face_neighbor.get_div_phi()));
   }
   if (!_unique_fe_face_neighbor_helper.empty())
   {
@@ -2981,6 +2999,8 @@ Assembly::copyShapes(unsigned int var)
     copyShapes(v);
     if (v.computingCurl())
       curlPhi(v).shallowCopy(v.curlPhi());
+    if (v.computingDiv())
+      divPhi(v).shallowCopy(v.divPhi());
   }
   else
     mooseError("Unsupported variable field type!");
@@ -3016,6 +3036,8 @@ Assembly::copyFaceShapes(unsigned int var)
     copyFaceShapes(v);
     if (v.computingCurl())
       _vector_curl_phi_face.shallowCopy(v.curlPhi());
+    if (v.computingDiv())
+      _vector_div_phi_face.shallowCopy(v.divPhi());
   }
   else
     mooseError("Unsupported variable field type!");
@@ -4709,6 +4731,42 @@ Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
   return _vector_fe_shape_data_face_neighbor[type]->_curl_phi;
 }
 
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhi<VectorValue<Real>>(FEType type) const
+{
+  _need_div[type] = true;
+  buildVectorFE(type);
+  return _vector_fe_shape_data[type]->_div_phi;
+}
+
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const
+{
+  _need_div[type] = true;
+  buildVectorFaceFE(type);
+  return _vector_fe_shape_data_face[type]->_div_phi;
+}
+
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const
+{
+  _need_div[type] = true;
+  buildVectorNeighborFE(type);
+  return _vector_fe_shape_data_neighbor[type]->_div_phi;
+}
+
+template <>
+const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
+Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
+{
+  _need_div[type] = true;
+  buildVectorFaceNeighborFE(type);
+  return _vector_fe_shape_data_face_neighbor[type]->_div_phi;
+}
+
 const MooseArray<ADReal> &
 Assembly::adCurvatures() const
 {
diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C
index b62e1e5c4c96..2f62326ba44c 100644
--- a/framework/src/variables/MooseVariableData.C
+++ b/framework/src/variables/MooseVariableData.C
@@ -47,6 +47,9 @@ MooseVariableData<OutputType>::MooseVariableData(const MooseVariableField<Output
     _need_curl(false),
     _need_curl_old(false),
     _need_curl_older(false),
+    _need_div(false),
+    _need_div_old(false),
+    _need_div_older(false),
     _need_ad(false),
     _need_ad_u(false),
     _need_ad_grad_u(false),
@@ -60,6 +63,8 @@ MooseVariableData<OutputType>::MooseVariableData(const MooseVariableField<Output
     _second_phi_face_assembly_method(nullptr),
     _curl_phi_assembly_method(nullptr),
     _curl_phi_face_assembly_method(nullptr),
+    _div_phi_assembly_method(nullptr),
+    _div_phi_face_assembly_method(nullptr),
     _ad_grad_phi_assembly_method(nullptr),
     _ad_grad_phi_face_assembly_method(nullptr),
     _time_integrator(nullptr),
@@ -92,6 +97,8 @@ MooseVariableData<OutputType>::MooseVariableData(const MooseVariableField<Output
       _second_phi_face_assembly_method = &Assembly::feSecondPhiFace<OutputShape>;
       _curl_phi_assembly_method = &Assembly::feCurlPhi<OutputShape>;
       _curl_phi_face_assembly_method = &Assembly::feCurlPhiFace<OutputShape>;
+      _div_phi_assembly_method = &Assembly::feDivPhi<OutputShape>;
+      _div_phi_face_assembly_method = &Assembly::feDivPhiFace<OutputShape>;
       _ad_grad_phi_assembly_method = &Assembly::feADGradPhi<OutputShape>;
       _ad_grad_phi_face_assembly_method = &Assembly::feADGradPhiFace<OutputShape>;
 
@@ -109,6 +116,8 @@ MooseVariableData<OutputType>::MooseVariableData(const MooseVariableField<Output
       _second_phi_face_assembly_method = &Assembly::feSecondPhiFaceNeighbor<OutputShape>;
       _curl_phi_assembly_method = &Assembly::feCurlPhiNeighbor<OutputShape>;
       _curl_phi_face_assembly_method = &Assembly::feCurlPhiFaceNeighbor<OutputShape>;
+      _div_phi_assembly_method = &Assembly::feDivPhiNeighbor<OutputShape>;
+      _div_phi_face_assembly_method = &Assembly::feDivPhiFaceNeighbor<OutputShape>;
 
       _ad_grad_phi = nullptr;
       _ad_grad_phi_face = nullptr;
@@ -155,6 +164,7 @@ MooseVariableData<OutputType>::setGeometry(Moose::GeometryType gm_type)
       _current_grad_phi = _grad_phi;
       _current_second_phi = _second_phi;
       _current_curl_phi = _curl_phi;
+      _current_div_phi = _div_phi;
       _current_ad_grad_phi = _ad_grad_phi;
       break;
     }
@@ -165,6 +175,7 @@ MooseVariableData<OutputType>::setGeometry(Moose::GeometryType gm_type)
       _current_grad_phi = _grad_phi_face;
       _current_second_phi = _second_phi_face;
       _current_curl_phi = _curl_phi_face;
+      _current_div_phi = _div_phi_face;
       _current_ad_grad_phi = _ad_grad_phi_face;
       break;
     }
@@ -329,6 +340,37 @@ MooseVariableData<OutputType>::curlSln(Moose::SolutionState state) const
   }
 }
 
+template <typename OutputType>
+const typename MooseVariableData<OutputType>::FieldVariableDivergence &
+MooseVariableData<OutputType>::divSln(Moose::SolutionState state) const
+{
+  divPhi();
+  divPhiFace();
+  switch (state)
+  {
+    case Moose::Current:
+    {
+      _need_div = true;
+      return _div_u;
+    }
+
+    case Moose::Old:
+    {
+      _need_div_old = true;
+      return _div_u_old;
+    }
+
+    case Moose::Older:
+    {
+      _need_div_older = true;
+      return _div_u_older;
+    }
+
+    default:
+      mooseError("We don't currently support divergence from the previous non-linear iteration");
+  }
+}
+
 template <typename OutputType>
 const typename MooseVariableData<OutputType>::FieldVariablePhiSecond &
 MooseVariableData<OutputType>::secondPhi() const
@@ -361,6 +403,22 @@ MooseVariableData<OutputType>::curlPhiFace() const
   return *_curl_phi_face;
 }
 
+template <typename OutputType>
+const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
+MooseVariableData<OutputType>::divPhi() const
+{
+  _div_phi = &_div_phi_assembly_method(_assembly, _fe_type);
+  return *_div_phi;
+}
+
+template <typename OutputType>
+const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
+MooseVariableData<OutputType>::divPhiFace() const
+{
+  _div_phi_face = &_div_phi_face_assembly_method(_assembly, _fe_type);
+  return *_div_phi_face;
+}
+
 template <typename OutputType>
 void
 MooseVariableData<OutputType>::computeValues()
@@ -392,6 +450,9 @@ MooseVariableData<OutputType>::computeValues()
   if (_need_curl)
     _curl_u.resize(nqp);
 
+  if (_need_div)
+    _div_u.resize(nqp);
+
   if (_need_second_previous_nl)
     _second_u_previous_nl.resize(nqp);
 
@@ -426,6 +487,9 @@ MooseVariableData<OutputType>::computeValues()
 
     if (_need_curl_old)
       _curl_u_old.resize(nqp);
+    
+    if (_need_div_old)
+      _div_u_old.resize(nqp);
 
     if (_need_second_older)
       _second_u_older.resize(nqp);
@@ -450,6 +514,9 @@ MooseVariableData<OutputType>::computeValues()
 
     if (_need_curl)
       _curl_u[i] = 0;
+    
+    if (_need_div)
+      _div_u[i] = 0;
 
     if (_need_second_previous_nl)
       _second_u_previous_nl[i] = 0;
@@ -488,12 +555,16 @@ MooseVariableData<OutputType>::computeValues()
 
       if (_need_curl_old)
         _curl_u_old[i] = 0;
+      
+      if (_need_div_old)
+        _div_u_old[i] = 0;
     }
   }
 
   bool second_required =
       _need_second || _need_second_old || _need_second_older || _need_second_previous_nl;
   bool curl_required = _need_curl || _need_curl_old;
+  bool div_required = _need_div || _need_div_old;
 
   for (unsigned int i = 0; i < num_dofs; i++)
   {
@@ -567,6 +638,19 @@ MooseVariableData<OutputType>::computeValues()
           _curl_u_old[qp] += curl_phi_local * _vector_tags_dof_u[_old_solution_tag][i];
       }
 
+      if (div_required)
+      {
+        mooseAssert(_current_div_phi,
+                    "We're requiring a divergence calculation but have not set a div shape function!");
+        const OutputShapeDivergence div_phi_local = (*_current_div_phi)[i][qp];
+
+        if (_need_div)
+          _div_u[qp] += div_phi_local * _vector_tags_dof_u[_solution_tag][i];
+
+        if (is_transient && _need_div_old)
+          _div_u_old[qp] += div_phi_local * _vector_tags_dof_u[_old_solution_tag][i];
+      }
+
       for (auto tag : _required_vector_tags)
       {
         if (_sys.hasVector(tag) && _sys.getVector(tag).closed())
@@ -646,6 +730,9 @@ MooseVariableData<RealEigenVector>::computeValues()
   if (_need_curl)
     _curl_u.resize(nqp);
 
+  if (_need_div)
+    _div_u.resize(nqp);
+
   if (_need_second_previous_nl)
     _second_u_previous_nl.resize(nqp);
 
@@ -681,6 +768,9 @@ MooseVariableData<RealEigenVector>::computeValues()
     if (_need_curl_old)
       _curl_u_old.resize(nqp);
 
+    if (_need_div_old)
+      _div_u_old.resize(nqp);
+
     if (_need_second_older)
       _second_u_older.resize(nqp);
   }
@@ -705,6 +795,9 @@ MooseVariableData<RealEigenVector>::computeValues()
     if (_need_curl)
       _curl_u[i].setZero(_count);
 
+    if (_need_div)
+      _div_u[i].setZero(_count);
+
     if (_need_second_previous_nl)
       _second_u_previous_nl[i].setZero(_count, LIBMESH_DIM * LIBMESH_DIM);
 
@@ -742,12 +835,16 @@ MooseVariableData<RealEigenVector>::computeValues()
 
       if (_need_curl_old)
         _curl_u_old[i].setZero(_count);
+      
+      if (_need_div_old)
+        _div_u_old[i].setZero(_count);
     }
   }
 
   bool second_required =
       _need_second || _need_second_old || _need_second_older || _need_second_previous_nl;
   bool curl_required = _need_curl || _need_curl_old;
+  bool div_required = _need_div || _need_div_old;
 
   for (unsigned int i = 0; i < num_dofs; i++)
   {
@@ -832,6 +929,19 @@ MooseVariableData<RealEigenVector>::computeValues()
           _curl_u_old[qp] += curl_phi_local * _vector_tags_dof_u[_old_solution_tag][i];
       }
 
+      if (div_required)
+      {
+        mooseAssert(_current_div_phi,
+                    "We're requiring a divergence calculation but have not set a divergence shape function!");
+        const auto div_phi_local = (*_current_div_phi)[i][qp];
+
+        if (_need_div)
+          _div_u[qp] += div_phi_local * _vector_tags_dof_u[_solution_tag][i];
+
+        if (is_transient && _need_div_old)
+          _div_u_old[qp] += div_phi_local * _vector_tags_dof_u[_old_solution_tag][i];
+      }
+
       for (auto tag : _required_vector_tags)
       {
         if (_need_vector_tag_u[tag])
diff --git a/framework/src/variables/MooseVariableFE.C b/framework/src/variables/MooseVariableFE.C
index 7c73bc13805f..c303f5915e8a 100644
--- a/framework/src/variables/MooseVariableFE.C
+++ b/framework/src/variables/MooseVariableFE.C
@@ -737,6 +737,13 @@ MooseVariableFE<OutputType>::curlPhi() const
   return _element_data->curlPhi();
 }
 
+template <typename OutputType>
+const typename MooseVariableFE<OutputType>::FieldVariablePhiDivergence &
+MooseVariableFE<OutputType>::divPhi() const
+{
+  return _element_data->divPhi();
+}
+
 template <typename OutputType>
 const typename MooseVariableFE<OutputType>::FieldVariablePhiSecond &
 MooseVariableFE<OutputType>::secondPhiFace() const
@@ -751,6 +758,13 @@ MooseVariableFE<OutputType>::curlPhiFace() const
   return _element_data->curlPhiFace();
 }
 
+template <typename OutputType>
+const typename MooseVariableFE<OutputType>::FieldVariablePhiDivergence &
+MooseVariableFE<OutputType>::divPhiFace() const
+{
+  return _element_data->divPhiFace();
+}
+
 template <typename OutputType>
 const typename MooseVariableFE<OutputType>::FieldVariablePhiSecond &
 MooseVariableFE<OutputType>::secondPhiNeighbor() const
@@ -765,6 +779,13 @@ MooseVariableFE<OutputType>::curlPhiNeighbor() const
   return _neighbor_data->curlPhi();
 }
 
+template <typename OutputType>
+const typename MooseVariableFE<OutputType>::FieldVariablePhiDivergence &
+MooseVariableFE<OutputType>::divPhiNeighbor() const
+{
+  return _neighbor_data->divPhi();
+}
+
 template <typename OutputType>
 const typename MooseVariableFE<OutputType>::FieldVariablePhiSecond &
 MooseVariableFE<OutputType>::secondPhiFaceNeighbor() const
@@ -779,6 +800,13 @@ MooseVariableFE<OutputType>::curlPhiFaceNeighbor() const
   return _neighbor_data->curlPhiFace();
 }
 
+template <typename OutputType>
+const typename MooseVariableFE<OutputType>::FieldVariablePhiDivergence &
+MooseVariableFE<OutputType>::divPhiFaceNeighbor() const
+{
+  return _neighbor_data->divPhiFace();
+}
+
 template <typename OutputType>
 bool
 MooseVariableFE<OutputType>::usesSecondPhi() const
@@ -800,6 +828,13 @@ MooseVariableFE<OutputType>::computingCurl() const
   return _element_data->computingCurl();
 }
 
+template <typename OutputType>
+bool
+MooseVariableFE<OutputType>::computingDiv() const
+{
+  return _element_data->computingDiv();
+}
+
 template <typename OutputType>
 bool
 MooseVariableFE<OutputType>::isNodalDefined() const