damage implemented in torch and fiber materials

src/mumfim/microscale/BatchedFiberRVEAnalysisExplicit.h

@@ -77,6 +77,7 @@ namespace mumfim
     ElementDataType current_length;
     // Only deal with fiber networks with uniform material properties
     PackedScalarType fiber_elastic_modulus;
+    PackedScalarType original_fiber_elastic_modulus;
     PackedScalarType fiber_area;
     PackedScalarType fiber_density;
     PackedScalarType viscous_damping_coefficient;
@@ -237,6 +238,8 @@ namespace mumfim
             fiber_elastic_modulus.template getRow<HostMemorySpace>(i);
         fiber_elastic_modulus_row(0) =
             fiber_networks[i]->getFiberReaction(0).getYoungModulus();
+        original_fiber_elastic_modulus.template getRow<HostMemorySpace>(i)(0) =
+            fiber_networks[i]->getFiberReaction(0).getYoungModulus();
         auto fiber_area_row = fiber_area.template getRow<HostMemorySpace>(i);
         fiber_area_row(0) =
             fiber_networks[i]->getFiberReaction(0).getFiberArea();
@@ -314,6 +317,7 @@ namespace mumfim
                                   current_deformation_gradient_,
                                   trial_deformation_gradient_);
       }
+
       // if we aren't "updating coords" we are doing a finite difference and
       // that should make use of the last state, not the accepted state!
       else
@@ -467,6 +471,14 @@ namespace mumfim
     auto GetStrainEnergy() {
       return strain_energy;
     }
+    void updateDamageFactor(Kokkos::View<Scalar*> damage_factor) override
+    {
+        assert(damage_factor.extent(0) == original_fiber_elastic_modulus.getNumRows());
+        Kokkos::parallel_for(damage_factor.extent(0), KOKKOS_LAMBDA(int i) {
+            auto original_modulus_d = original_fiber_elastic_modulus.template getRow<DeviceMemorySpace>(i);
+            fiber_elastic_modulus.template getRow<DeviceMemorySpace>(i)(0) = original_modulus_d(0) * damage_factor(i);
+        });
+    }
 
     private:
     BatchedAnalysisGetPK2StressFunc<BatchedFiberRVEAnalysisExplicit>

src/mumfim/microscale/BatchedRVEAnalysis.h

@@ -56,6 +56,9 @@ namespace mumfim
       Kokkos::deep_copy(omega.h_view, 0);
     }
     size_t GetNumRVEs() const { return this->num_rves_; }
+
+    virtual void updateDamageFactor(Kokkos::View<Scalar*> damage_factor) {}
+
     // RVEAnalysis(const RVEAnalysis & an);
     // RVEAnalysis();
   };

src/mumfim/microscale/BatchedTorchAnalysis.h

@@ -6,6 +6,7 @@
 #include <mumfim/microscale/ContinuumMechanics.h>
 #include <mumfim/torch/TorchUtilities.h>
 #include <torch/script.h>
+#include <KokkosBlas1_scal.hpp>
 
 #include <Kokkos_Core.hpp>
 
@@ -25,9 +26,11 @@ namespace mumfim
     BatchedTorchAnalysis(const std::string & model_path, int nrves)
         : BatchedRVEAnalysis<Scalar, LocalOrdinal, ExeSpace>(nrves),
           F_trial_("F_trial", nrves), stiffness_trial_("stiffness_trial", nrves),
-          F_current_("F_current", nrves)
+          F_current_("F_current", nrves),
+          damage_factor_("damage_factor", nrves)
     {
       FillIdentity(F_current_);
+      Kokkos::deep_copy(damage_factor_, 1.0);
       try
       {
         model_ = torch::jit::load(model_path);
@@ -66,6 +69,16 @@ namespace mumfim
                                            right_cauchy_green_[0].toTensor());
       ml::TorchArrayToKokkosView(cauchy_stress, sigma.template view<ExeSpace>());
       ml::TorchArrayToKokkosView(stiffness, stiffness_trial_);
+
+      // multiply stress an stiffness by damage
+      //KokkosBlas::scal(stiffness_trial_, damage_factor_, stiffness_trial_);
+      KokkosBlas::scal(sigma.d_view, damage_factor_, sigma.d_view);
+      // since stiffness is 3D, need to do team loop...
+      Kokkos::MDRangePolicy policy({0,0,0},{stiffness_trial_.extent(0),stiffness_trial_.extent(1),stiffness_trial_.extent(2)});
+      Kokkos::parallel_for(policy,
+                           KOKKOS_LAMBDA(const int i, const int j, const int k) {
+                             stiffness_trial_(i, j,k) *= damage_factor_(i);
+                           });
       sigma.template modify<ExeSpace>();
 
       return true;
@@ -85,6 +98,10 @@ namespace mumfim
       KOKKOS_ASSERT(C.extent(0) == this->num_rves_);
       Kokkos::deep_copy(C.template view<ExeSpace>(), stiffness_trial_);
     }
+    void updateDamageFactor(Kokkos::View<Scalar*> damage_factor) override
+    {
+      Kokkos::deep_copy(damage_factor_, damage_factor);
+    }
 
     private:
     using LO = LocalOrdinal;
@@ -96,6 +113,7 @@ namespace mumfim
 
     // current deformation gradient that is updated after each accept
     Kokkos::View<Scalar * [3][3], ExeSpace> F_current_;
+    Kokkos::View<Scalar* > damage_factor_;
 
     torch::jit::script::Module model_;
 

src/mumfim/microscale/MultiscaleRVEAnalysis.cc

@@ -25,6 +25,21 @@
 #endif
 namespace mumfim
 {
+
+  static void calculateUpdatedDamageFactor(double failure_stress, double damage_factor,
+    const Kokkos::View<Scalar * [6]> & stress, const Kokkos::View<Scalar*> & accepted_damage,
+    const Kokkos::View<Scalar*> & trial_damage)
+  {
+    auto von_mises = computeVonMisesStress(stress);
+    Kokkos::parallel_for("Calculate Updated Damage Factor",
+      Kokkos::RangePolicy<>(0, trial_damage.size()),
+      KOKKOS_LAMBDA(const int i) {
+            trial_damage(i) = von_mises(i) > failure_stress
+                                  ? accepted_damage(i) * damage_factor
+                                  : accepted_damage(i);
+      });
+  }
+
   MultiscaleRVEAnalysis::~MultiscaleRVEAnalysis() = default;
 
   MultiscaleRVEAnalysis::MultiscaleRVEAnalysis(
@@ -292,9 +307,10 @@ namespace mumfim
                                                          hdrs.size());
     Kokkos::DualView<Scalar * [3]> orientation_tensor_normal(
         "orientation tensor normal", hdrs.size());
-    Kokkos::View<Scalar *, Kokkos::HostSpace> trial_damage("damage factor",
+    Kokkos::View<Scalar *> trial_damage("damage factor",
                                                            hdrs.size());
-    Kokkos::View<Scalar *, Kokkos::HostSpace> accepted_damage("damage factor",
+    auto trial_damage_h = Kokkos::create_mirror_view(trial_damage);
+    Kokkos::View<Scalar *> accepted_damage("damage factor",
                                                               hdrs.size());
     Kokkos::deep_copy(accepted_damage, 1.0);
     // communicate the step results back to the macro scale
@@ -313,7 +329,6 @@ namespace mumfim
       int step_accepted{0};
       while (!step_complete)
       {
-        Kokkos::deep_copy(trial_damage, accepted_damage);
         // migration
         // if (macro_iter == 0) updateCoupling();
         // send the initial microscale rve states back to the macroscale
@@ -329,10 +344,6 @@ namespace mumfim
         auto stress_h = stress.h_view;
         auto material_stiffness_h = material_stiffness.h_view;
         deformation_gradient.modify<Kokkos::HostSpace>();
-        for (size_t i = 0; i < trial_damage.size(); ++i)
-        {
-          trial_damage(i) *= damage_factor_;
-        }
         // fill the deformation gradient data
         for (std::size_t i = 0; i < results.size(); ++i)
         {
@@ -346,15 +357,19 @@ namespace mumfim
         }
         batched_analysis->run(deformation_gradient, stress);
         batched_analysis->computeMaterialStiffness(material_stiffness);
+        calculateUpdatedDamageFactor(failure_stress_, damage_factor_, stress.d_view, accepted_damage,
+                                     trial_damage);
+        batched_analysis->updateDamageFactor(trial_damage);
         stress.sync<Kokkos::HostSpace>();
         material_stiffness.sync<Kokkos::HostSpace>();
+        Kokkos::deep_copy(trial_damage_h, trial_damage);
         // fill the results data
         for (std::size_t i = 0; i < results.size(); ++i)
         {
           double * sigma = &(results[i].data[0]);
           double * avgVolStress = &(results[i].data[6]);
           double * matStiffMatrix = &(results[i].data[9]);
-          results[i].data[45] = trial_damage(i);
+          results[i].data[45] = trial_damage_h(i);
           for (int j = 0; j < 6; ++j)
           {
             sigma[j] = stress_h(i, j);

src/mumfim/microscale/batched_analysis/BatchedFiberRVEAnalysisExplicitBase.h

@@ -293,5 +293,18 @@ namespace mumfim
       }
     }
   };
+  inline Kokkos::View<Scalar*> computeVonMisesStress(Kokkos::View<Scalar*[6]> stress) {
+    Kokkos::View<Scalar*> vonMisesStress("von mises stress", stress.extent(0));
+    Kokkos::parallel_for(stress.extent(0), KOKKOS_LAMBDA(int i) { 
+        Scalar sigma_v_sq = 0.5*( pow(stress(i,0)-stress(i,1),2) +
+                                  pow(stress(i,1)-stress(i,2),2) +
+                                  pow(stress(i,2)-stress(i,0),2) +
+                                  6*(pow(stress(i,3),2) +
+                                     pow(stress(i,4),2) +
+                                     pow(stress(i,5),2)));
+        vonMisesStress(i) = sqrt(sigma_v_sq);
+        });
+    return vonMisesStress;
+  }
 }  // namespace mumfim
 #endif