Feature/fsisphInterfaces

RELEASE_NOTES.md

@@ -19,7 +19,12 @@ Notable changes include:
   * Bug Fixes / improvements:
     * Fixed melt behavior in Steinberg-Guinan strength model, which was ignoring melt for damaged material.
     * Fixed range of dimensionless melt temperature for Johnson-Cook strength.
-
+    * FSISPH new features and modifications to method 
+      * NOTE constructor inputs have changed
+      * strength implementation modified
+      * new features added including plane strain option and settable minP for interfaces
+      * new, more rigorous, interface and free surface tracking
+
 Version v2023-06-0 -- Release date 2023-06-20
 ==============================================
 

src/DEM/DEMBase.cc

@@ -19,6 +19,7 @@
 #include "DataBase/ReplaceFieldList.hh"
 #include "DataBase/DataBase.hh"
 
+
 #include "Field/FieldList.hh"
 #include "Field/NodeIterators.hh"
 

src/Damage/DamagedPressurePolicy.cc

@@ -4,6 +4,7 @@
 //----------------------------------------------------------------------------//
 #include "DamagedPressurePolicy.hh"
 #include "Hydro/HydroFieldNames.hh"
+#include "FSISPH/FSIFieldNames.hh"
 #include "Strength/SolidFieldNames.hh"
 #include "NodeList/SolidNodeList.hh"
 #include "SolidMaterial/SolidEquationOfState.hh"
@@ -51,7 +52,8 @@ update(const KeyType& key,
        const double dt) {
   KeyType fieldKey, nodeListKey, Dkey;
   StateBase<Dimension>::splitFieldKey(key, fieldKey, nodeListKey);
-  REQUIRE(fieldKey == HydroFieldNames::pressure and 
+  REQUIRE((fieldKey == HydroFieldNames::pressure or
+           fieldKey == FSIFieldNames::damagedPressure) and 
           nodeListKey == UpdatePolicyBase<Dimension>::wildcard());
   StateBase<Dimension>::buildFieldKey(SolidFieldNames::tensorDamage, nodeListKey);
   auto pressure = state.fields(fieldKey, 0.0);

src/Damage/TensorStrainPolicy.cc

@@ -163,8 +163,20 @@ update(const KeyType& key,
       }
     }
 
+    //------------------------------------------------------------------------------------
+    // NOTE: this is a temporary fix to address the difference between Spheral's damaged
+    //       pressure and the standard EOS pressure originally used by Benz and Asphaug.
+    //       (1-D) is squared to counter-act the 1-D rolled into the pressure. Original 
+    //        code is commented out below. (This will enhance the deviatoric portion)
+    //       -JMPearl 11/14/2023
+
+    //const auto fDs = 1.0 - D(i).eigenValues().maxElement();
+    //stateField(i) *= safeInvVar(max(0.0, fDs*fDs), tiny);
+
     // Damage enhancement of the effective strain.
     stateField(i) *= safeInvVar(max(0.0, 1.0 - D(i).Trace()/Dimension::nDim), tiny);
+    //------------------------------------------------------------------------------------
+
 
     // Apply limiting to the effective strain.
     stateField(i) = max(1.0e-7*max(1.0, std::abs(stateField(i).Trace())/Dimension::nDim), stateField(i));

src/FSISPH/CMakeLists.txt

@@ -1,8 +1,10 @@
 include_directories(.)
+
 set(FSISPH_inst
     SolidFSISPHHydroBase
     computeFSISPHSumMassDensity
     computeHWeightedFSISPHSumMassDensity
+    computeInterfacePressureCorrectedSumMassDensity
     SlideSurface
     )
 
@@ -18,6 +20,7 @@ set(FSISPH_headers
     SolidFSISPHHydroBaseInline.hh
     computeFSISPHSumMassDensity.hh
     computeHWeightedFSISPHSumMassDensity.hh
+    computeInterfacePressureCorrectedSumMassDensity.hh
     FSIFieldNames.hh
     SlideSurface.hh
     SlideSurfaceInline.hh

src/FSISPH/FSIFieldNames.cc

@@ -6,11 +6,14 @@
 
 #include "FSIFieldNames.hh"
 
-const std::string Spheral::FSIFieldNames::rawPressure = "raw eos pressure";
+const std::string Spheral::FSIFieldNames::damagedPressure = "damaged pressure";
 const std::string Spheral::FSIFieldNames::pressureGradient = "pressureGradient";
 const std::string Spheral::FSIFieldNames::specificThermalEnergyGradient = "specificThermalEnergyGradient";
+const std::string Spheral::FSIFieldNames::interfaceFlags = "interfaceFlags";
+const std::string Spheral::FSIFieldNames::interfaceAreaVectors = "interfaceAreaVectors";
 const std::string Spheral::FSIFieldNames::interfaceNormals = "interfaceNormals";
+const std::string Spheral::FSIFieldNames::interfaceAngles = "interfaceAngles";
 const std::string Spheral::FSIFieldNames::interfaceFraction = "interfaceFraction";
 const std::string Spheral::FSIFieldNames::interfaceSmoothness = "interfaceSmoothness";
 const std::string Spheral::FSIFieldNames::smoothedInterfaceNormals = "smoothedInterfaceNormals";
-const std::string Spheral::FSIFieldNames::smoothnessNormalization = "smoothnessNormalization";
+const std::string Spheral::FSIFieldNames::interfaceSmoothnessNormalization = "interfaceSmoothnessNormalization";

src/FSISPH/FSIFieldNames.hh

@@ -9,14 +9,18 @@
 namespace Spheral {
 
 struct FSIFieldNames {
-  static const std::string rawPressure;
+  static const std::string damagedPressure;
   static const std::string pressureGradient;
   static const std::string specificThermalEnergyGradient;
+  static const std::string interfaceFlags;
+  static const std::string interfaceAreaVectors;
   static const std::string interfaceNormals;
+  static const std::string interfaceAngles;
   static const std::string interfaceFraction;
   static const std::string interfaceSmoothness;
   static const std::string smoothedInterfaceNormals;
-  static const std::string smoothnessNormalization;
+  static const std::string interfaceSmoothnessNormalization;
+  //static const std::string inverseEquivalentDeviatoricStress;
 };
 
 }

src/FSISPH/FSISPHHydros.py

@@ -7,7 +7,6 @@ def FSISPH(dataBase,
         W,
         Q = None,
         slides=None,
-        filter = 0.0,
         cfl = 0.35,
         surfaceForceCoefficient=0.0,
         densityStabilizationCoefficient=0.1, 
@@ -19,12 +18,14 @@ def FSISPH(dataBase,
         useVelocityMagnitudeForDt = False,
         compatibleEnergyEvolution = True,
         evolveTotalEnergy = False,
-        correctVelocityGradient = True,    
+        linearCorrectGradients = True,
+        planeStrain = False,
+        interfacePmin = 0.0,
+        interfaceNeighborAngleThreshold=0.707,
         HUpdate = IdealH,
+        densityUpdate = FSISumMassDensity,
         epsTensile = 0.0,
         nTensile = 4.0,
-        damageRelieveRubble = False,
-        strengthInDamage = False,
         xmin = (-1e100, -1e100, -1e100),
         xmax = ( 1e100,  1e100,  1e100),
         ASPH = False,
@@ -38,9 +39,6 @@ def FSISPH(dataBase,
     if compatibleEnergyEvolution and evolveTotalEnergy:
         raise RuntimeError("compatibleEnergyEvolution and evolveTotalEnergy are incompatible")
 
-    if strengthInDamage and damageRelieveRubble:
-        raise RuntimeError("strengthInDamage and damageRelieveRubble are incompatible")
-
     # create the map nodelist --> index
     nodeLists = dataBase.nodeLists()
     nodeListMap = {}
@@ -101,7 +99,6 @@ def FSISPH(dataBase,
               "Q" : Q,
               "slides" : slides,
               "W" : W,
-              "filter" : filter,
               "cfl" : cfl,
               "surfaceForceCoefficient" : surfaceForceCoefficient,
               "densityStabilizationCoefficient" : densityStabilizationCoefficient,      
@@ -113,15 +110,14 @@ def FSISPH(dataBase,
               "useVelocityMagnitudeForDt" : useVelocityMagnitudeForDt,
               "compatibleEnergyEvolution" : compatibleEnergyEvolution,
               "evolveTotalEnergy" : evolveTotalEnergy,
-              "gradhCorrection" : False,
-              "XSPH" : False,
-              "correctVelocityGradient" : correctVelocityGradient,
-              "densityUpdate" : IntegrateDensity,
+              "linearCorrectGradients" : linearCorrectGradients,
+              "planeStrain" : planeStrain,
+              "interfacePmin" : interfacePmin,
+              "interfaceNeighborAngleThreshold" : interfaceNeighborAngleThreshold,
+              "densityUpdate" : densityUpdate,
               "HUpdate" : HUpdate,
               "epsTensile" : epsTensile,
               "nTensile" : nTensile,
-              "damageRelieveRubble" : damageRelieveRubble,
-              "strengthInDamage" : strengthInDamage,
               "xmin" : eval("Vector%id(%g, %g, %g)" % xmin),
               "xmax" : eval("Vector%id(%g, %g, %g)" % xmax)}