Merge pull request #317 from LLNL/feature/ArtificialViscosityPackage

Splitting Artificial viscosity out of hydrodynamics packages into its own physics package

Author: jmikeowen

RELEASE_NOTES.md

@@ -34,6 +34,11 @@ Notable changes include:
       is correct.
     * Performance regression testing is now available. All developers are encouraged to run the performance testing suite for any code changes that might impact performance. See documentation for more details.
     * Added our old ASPH IdealH H update as an option. While it is not as reliable as our current default ASPH, it does not require building the Voronoi and is therefore signifcantly faster.
+    * Converted artificial viscosities to Physics packages, and add them as pre-subpackages to Hydro objects.
+    * Split artificial viscosities based on the type of pressure they compute (currently Scalar or Tensor), which is slightly more efficient.
+      * This required making the hydro packages evaluateDerivatives into templated methods based on the type of Q they are handed.
+      * Also introduced a new base class (ArtificialViscosityHandle), which provides a handle class not templated on the type
+        of Q pressure for Hydro objects to hold onto.
 
   * Build changes / improvements:
     * Distributed source directory must always be built now.

scripts/devtools/performance_analysis.py

@@ -255,7 +255,7 @@ def get_caliper_files(file_path):
         th.stats.std(rtest, metrics)
     # Compute the max allowable time for the main region
     ctest.statsframe.dataframe["thresh"] = compute_threshold(rtest.statsframe, metric0)
-    ctest.move_metrics_to_statsframe([metric0])
+    ctest.move_metrics_to_statsframe([metric0,metric1])
     ref_main = get_times(rtest.statsframe, "main", metric0+"_mean")[0]
     cur_main = get_times(ctest.statsframe, "main", metric0)[0]
     ref_thresh = get_times(ctest.statsframe, "main", "thresh")[0]

src/ArtificialViscosity/ArtificialViscosity.cc

@@ -2,230 +2,74 @@
 // ArtificialViscosity -- The base class for all ArtificialViscosities in 
 // Spheral++.
 //
+// This intermediate class between ArtificialViscosityHandle and the actual
+// artficial viscosity implementations specifies the QPi = Q/rho^2 return type,
+// generally either a Scalar or a Tensor.
+//
 // Created by JMO, Sun May 21 21:16:43 PDT 2000
 //----------------------------------------------------------------------------//
 #ifndef __Spheral_ArtificialViscosity__
 #define __Spheral_ArtificialViscosity__
 
-#include <utility>
-
-#include "Geometry/Dimension.hh"
-#include "RK/RKCorrectionParams.hh"
-
-#include <vector>
-#include "Field/FieldList.hh"
-#include "DataOutput/registerWithRestart.hh"
-
-namespace Spheral {
-
-  template<typename Dimension> class AllNodeIterator;
-  template<typename Dimension> class InternalNodeIterator;
-  template<typename Dimension> class GhostNodeIterator;
-  template<typename Dimension> class MasterNodeIterator;
-  template<typename Dimension> class CoarseNodeIterator;
-  template<typename Dimension> class RefineNodeIterator;
-
-  template<typename Dimension> class State;
-  template<typename Dimension> class StateDerivatives;
+#include "ArtificialViscosity/ArtificialViscosityHandle.hh"
 
-  template<typename Dimension> class DataBase;
-  template<typename Dimension, typename DataType> class FieldList;
-  class FileIO;
-  template<typename Dimension> class ConnectivityMap;
-  template<typename Dimension> class TableKernel;
-  template<typename Dimension> class Boundary;
-}
+#include <utility>
+#include <typeindex>
 
 namespace Spheral {
 
-template<typename Dimension>
-class ArtificialViscosity {
+template<typename Dimension, typename QPiType>
+class ArtificialViscosity: public ArtificialViscosityHandle<Dimension> {
 public:
   //--------------------------- Public Interface ---------------------------//
-  typedef typename Dimension::Scalar Scalar;
-  typedef typename Dimension::Vector Vector;
-  typedef typename Dimension::Tensor Tensor;
-  typedef typename Dimension::SymTensor SymTensor;
-  typedef typename std::pair<double, std::string> TimeStepType;
+  using Scalar = typename Dimension::Scalar;
+  using Vector = typename Dimension::Vector;
+  using Tensor = typename Dimension::Tensor;
+  using SymTensor = typename Dimension::SymTensor;
 
-  typedef typename std::vector<Boundary<Dimension>*>::const_iterator ConstBoundaryIterator;
+  using ReturnType = QPiType;
 
-  // Constructors.
+  // Constructors, destructor
   ArtificialViscosity(const Scalar Clinear,
                       const Scalar Cquadratic,
-                      const RKOrder QcorrectionOrder = RKOrder::LinearOrder);
-
-  // Destructor.
-  virtual ~ArtificialViscosity();
-
-  // Initialize the artificial viscosity for all FluidNodeLists in the given
-  // DataBase.
-  virtual void initialize(const DataBase<Dimension>& dataBase,
-                          const State<Dimension>& state,
-                          const StateDerivatives<Dimension>& derivs,
-                          ConstBoundaryIterator boundaryBegin,
-                          ConstBoundaryIterator boundaryEnd,
-                          const Scalar /*time*/,
-                          const Scalar /*dt*/,
-                          const TableKernel<Dimension>& W);
-
-  // Require all descendents to return the artificial viscous Pi = P/rho^2 as a tensor.
-  // Scalar viscosities should just return a diagonal tensor with their value along the diagonal.
-  virtual std::pair<Tensor, Tensor> Piij(const unsigned nodeListi, const unsigned i, 
-                                         const unsigned nodeListj, const unsigned j,
-                                         const Vector& xi,
-                                         const Vector& etai,
-                                         const Vector& vi,
-                                         const Scalar rhoi,
-                                         const Scalar csi,
-                                         const SymTensor& Hi,
-                                         const Vector& xj,
-                                         const Vector& etaj,
-                                         const Vector& vj,
-                                         const Scalar rhoj,
-                                         const Scalar csj,
-                                         const SymTensor& Hj) const = 0;
-
-  // Allow access to the linear and quadratic constants.
-  Scalar Cl() const;
-  void Cl(Scalar Cl);
-
-  Scalar Cq() const;
-  void Cq(Scalar Cq);
-
-  //Allow access to the Q correction order.
-  RKOrder QcorrectionOrder() const;
-  void QcorrectionOrder(RKOrder order);
-
-  // Toggle for the Balsara shear correction.
-  bool balsaraShearCorrection() const;
-  void balsaraShearCorrection(bool value);
-
-  // Calculate the curl of the velocity given the stress tensor.
-  Scalar curlVelocityMagnitude(const Tensor& DvDx) const;
-
-  // Access the FieldList of multiplicative corrections.
-  FieldList<Dimension, Scalar>& ClMultiplier();
-  FieldList<Dimension, Scalar>& CqMultiplier();
-  FieldList<Dimension, Scalar>& shearCorrection();
-  const FieldList<Dimension, Scalar>& ClMultiplier() const;
-  const FieldList<Dimension, Scalar>& CqMultiplier() const;
-  const FieldList<Dimension, Scalar>& shearCorrection() const;
-
-  // Access the internally computed estimate of sigma:
-  // sig^ab = \partial v^a / \partial x^b.
-  const FieldList<Dimension, Tensor>& sigma() const;
-
-  // Access the internally computed estimate of the velocity gradient and
-  // grad div velocity.
-  const FieldList<Dimension, Vector>& gradDivVelocity() const;
-
-  // Switch to turn the del^2 v limiter on/off.
-  bool limiter() const;
-  void limiter(bool value);
-
-  // Access the epsilon safety factor.
-  Scalar epsilon2() const;
-  void epsilon2(Scalar epsilon2);
-
-  // Method to return the limiter magnitude for the given node.
-  Tensor calculateLimiter(const Vector& /*vi*/,
-                          const Vector& /*vj*/,
-                          const Scalar  ci,
-                          const Scalar  /*cj*/,
-                          const Scalar  hi,
-                          const Scalar  /*hj*/,
-                          const int nodeListID,
-                          const int nodeID) const;
-
-  // Helper for the limiter, calculate the unit grad div v term for the given 
-  // node.
-  Vector shockDirection(const Scalar ci,
-                        const Scalar hi,
-                        const int nodeListID,
-                        const int nodeID) const;
-
-  // The negligible sound speed parameter for use in the limiter.
-  Scalar negligibleSoundSpeed() const;
-  void negligibleSoundSpeed(Scalar val);
-
-  // The multiplier for sound speed in the limiter.
-  Scalar csMultiplier() const;
-  void csMultiplier(Scalar val);
-
-  // The multiplier for energy in the limiter.
-  Scalar energyMultiplier() const;
-  void energyMultiplier(Scalar val);
-
-  // Helper method to calculate Del cross V from the given sigma tensor.
-  Scalar computeDelCrossVMagnitude(const Tensor& sigma) const;
-
-  // Helper method to calculate the weighting based on the given position
-  // for use in the sigma calculation.
-  Vector sigmaWeighting(const Vector& r) const;
-
-  // Figure out the total stress-strain tensor for a given node pair based on 
-  // the stored average value and the given (position, velocity) pair.
-  Tensor sigmaij(const Vector& rji, 
-                 const Vector& rjiUnit, 
-                 const Vector& vji, 
-                 const Scalar& hi2,
-                 const int nodeListID,
-                 const int nodeID) const;
-
-  //****************************************************************************
-  // Methods required for restarting.
-  virtual std::string label() const { return "ArtificialViscosity"; }
-  virtual void dumpState(FileIO& file, const std::string& pathName) const;
-  virtual void restoreState(const FileIO& file, const std::string& pathName);
-  //****************************************************************************
-
-  // Allow descendents to request that sigma be calculated.
-  bool calculateSigma() const;
-  void calculateSigma(bool value);
-
-protected:
-  //--------------------------- Protected Interface ---------------------------//
-  Scalar mClinear;
-  Scalar mCquadratic;
-  RKOrder mQcorrectionOrder;
-
-  // Switch for the Balsara shear correction.
-  bool mBalsaraShearCorrection;
-
-  // Generic multipliers for the linear and quadratic terms.
-  FieldList<Dimension, Scalar> mClMultiplier, mCqMultiplier, mShearCorrection;
-
-  // Parameters for the Q limiter.
-  bool mCalculateSigma;
-  bool mLimiterSwitch;
-  Scalar mEpsilon2;
-  Scalar mNegligibleSoundSpeed;
-  Scalar mCsMultiplier;
-  Scalar mEnergyMultiplier;
-  FieldList<Dimension, Tensor> mSigma;
-  FieldList<Dimension, Vector> mGradDivVelocity;
-    
-  // The restart registration.
-  RestartRegistrationType mRestart;
-
-  // Protected methods.
-  virtual void calculateSigmaAndGradDivV(const DataBase<Dimension>& dataBase,
-                                         const State<Dimension>& state,
-                                         const StateDerivatives<Dimension>& /*derivs*/,
-                                         const TableKernel<Dimension>& W,
-                                         ConstBoundaryIterator boundaryBegin,
-                                         ConstBoundaryIterator boundaryEnd);
-
-private:
-  //--------------------------- Private Interface ---------------------------//
-  ArtificialViscosity();
-  ArtificialViscosity(const ArtificialViscosity&);
-  ArtificialViscosity& operator=(const ArtificialViscosity&) const;
+                      const TableKernel<Dimension>& kernel): ArtificialViscosityHandle<Dimension>(Clinear, Cquadratic, kernel) {}
+  virtual ~ArtificialViscosity() = default;
+
+  // No default constructor, copying, or assignment
+  ArtificialViscosity() = delete;
+  ArtificialViscosity(const ArtificialViscosity&) = delete;
+  ArtificialViscosity& operator=(const ArtificialViscosity&) = delete;
+
+  //...........................................................................
+  // Virtual methods we expect ArtificialViscosities to provide
+  // Required method returning the type_index of the descendant QPiType
+  virtual std::type_index QPiTypeIndex() const override { return std::type_index(typeid(QPiType)); }
+
+  // All ArtificialViscosities must provide the pairwise QPi term (pressure/rho^2)
+  // Returns the pair values QPiij and QPiji by reference as the first two arguments.
+  // Note the final FieldLists (fCl, fCQ, DvDx) should be the special versions registered
+  // by the ArtficialViscosity (particularly DvDx).
+  virtual void QPiij(QPiType& QPiij, QPiType& QPiji,    // result for QPi (Q/rho^2)
+                     Scalar& Qij, Scalar& Qji,          // result for viscous pressure
+                     const unsigned nodeListi, const unsigned i, 
+                     const unsigned nodeListj, const unsigned j,
+                     const Vector& xi,
+                     const SymTensor& Hi,
+                     const Vector& etai,
+                     const Vector& vi,
+                     const Scalar rhoi,
+                     const Scalar csi,
+                     const Vector& xj,
+                     const SymTensor& Hj,
+                     const Vector& etaj,
+                     const Vector& vj,
+                     const Scalar rhoj,
+                     const Scalar csj,
+                     const FieldList<Dimension, Scalar>& fCl,
+                     const FieldList<Dimension, Scalar>& fCq,
+                     const FieldList<Dimension, Tensor>& DvDx) const = 0;
 };
 
 }
 
-#include "ArtificialViscosityInline.hh"
-
 #endif