From ff76b83c989d1d5ac850b36b9c5c09d9aa19ef6a Mon Sep 17 00:00:00 2001
From: timovanopstal <timovanopstal@gmail.com>
Date: Mon, 26 Sep 2016 18:52:09 +0200
Subject: [PATCH] add: variations of stabilization params w.r.t. velocity

---
 Apps/Common/StabilizationUtils.C | 19 +++++++++++++++++++
 Apps/Common/StabilizationUtils.h | 19 ++++++++++++++++++-
 2 files changed, 37 insertions(+), 1 deletion(-)

diff --git a/Apps/Common/StabilizationUtils.C b/Apps/Common/StabilizationUtils.C
index 4daf4247..e193bf48 100644
--- a/Apps/Common/StabilizationUtils.C
+++ b/Apps/Common/StabilizationUtils.C
@@ -59,4 +59,23 @@ bool getTauNSALEPt(double dt, double mu, const Vector& U, const Matrix& G, doubl
   return true;
 }
 
+
+bool getTauPtJac(const Vector& U, const Matrix& G, const double tauM, Vector& tauMjac)
+{
+  tauMjac = -pow(tauM,3)*G*U;
+
+  return true;
+}
+
+
+bool getTauNSPtJac(const Vector& U, const Matrix& G, const double tauM, const double& tauC, Vector& tauMjac, Vector& tauCjac)
+{
+  if (!getTauPtJac(U, G, tauM, tauMjac))
+    return false;
+
+  tauCjac = (-tauC/tauM) * tauMjac;
+
+  return true;
+}
+
 }
diff --git a/Apps/Common/StabilizationUtils.h b/Apps/Common/StabilizationUtils.h
index ac3dbd63..8350952b 100644
--- a/Apps/Common/StabilizationUtils.h
+++ b/Apps/Common/StabilizationUtils.h
@@ -21,7 +21,7 @@ namespace StabilizationUtils {
   //! \details Stabilization parameter in integration point
   double getTauPt(double dt, double mu, const Vector& U, const Matrix& G,
       const double Ct=2.0, const double Cl=36.0);
- 		  
+
   //! \brief Computes stabilization parameters for Navier-Stokes equations
   //! \param[in] dt The timestep size
   //! \param[in] mu Diffusion/viscosity parameter
@@ -45,6 +45,23 @@ namespace StabilizationUtils {
   //! \details Stabilization parameters in integration point
   bool getTauNSALEPt(double dt, double mu, const Vector& U, const Matrix& G,
 		     double& tauM, double& tauC, const double Ct=2.0, const double Cl=36.0);
+
+  //! \brief Computes variation of stability parameters with respect to velocity
+  //! \param[in] U  Velocity vector
+  //! \param[in] G  The G matrix
+  //! \param[in] tauM Stabilization parameter for momentum
+  //! \param[out] tauMjac Variation of tauM with respect to U
+  bool getTauPtJac(const Vector& U, const Matrix& G, const double tauM, Vector& tauMjac);
+
+  //! \brief Computes variation of stability parameters with respect to velocity
+  //! \param[in] U  Velocity vector
+  //! \param[in] G  The G matrix
+  //! \param[in] tauM Stabilization parameter for momentum
+  //! \param[in] tauC Stabilization parameter for continuity
+  //! \param[out] tauMjac Variation of tauM with respect to U
+  //! \param[out] tauCjac Variation of tauC with respect to U
+  bool getTauNSPtJac(const Vector& U, const Matrix& G, const double tauM,
+      const double& tauC, Vector& tauMjac, Vector& tauCjac);
 }
 
 #endif