From 4580ffc972b289887a07bf5dc7dcd8d37a91f0b1 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 27 Jul 2024 16:38:22 -0700 Subject: [PATCH] damage implemented in torch and fiber materials --- .../BatchedFiberRVEAnalysisExplicit.h | 12 +++++++ src/mumfim/microscale/BatchedRVEAnalysis.h | 3 ++ src/mumfim/microscale/BatchedTorchAnalysis.h | 20 +++++++++++- .../microscale/MultiscaleRVEAnalysis.cc | 31 ++++++++++++++----- .../BatchedFiberRVEAnalysisExplicitBase.h | 13 ++++++++ 5 files changed, 70 insertions(+), 9 deletions(-) diff --git a/src/mumfim/microscale/BatchedFiberRVEAnalysisExplicit.h b/src/mumfim/microscale/BatchedFiberRVEAnalysisExplicit.h index f0c12b9f..95539bd5 100644 --- a/src/mumfim/microscale/BatchedFiberRVEAnalysisExplicit.h +++ b/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(i); fiber_elastic_modulus_row(0) = fiber_networks[i]->getFiberReaction(0).getYoungModulus(); + original_fiber_elastic_modulus.template getRow(i)(0) = + fiber_networks[i]->getFiberReaction(0).getYoungModulus(); auto fiber_area_row = fiber_area.template getRow(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 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(i); + fiber_elastic_modulus.template getRow(i)(0) = original_modulus_d(0) * damage_factor(i); + }); + } private: BatchedAnalysisGetPK2StressFunc diff --git a/src/mumfim/microscale/BatchedRVEAnalysis.h b/src/mumfim/microscale/BatchedRVEAnalysis.h index 108c4dbc..a756cf32 100644 --- a/src/mumfim/microscale/BatchedRVEAnalysis.h +++ b/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 damage_factor) {} + // RVEAnalysis(const RVEAnalysis & an); // RVEAnalysis(); }; diff --git a/src/mumfim/microscale/BatchedTorchAnalysis.h b/src/mumfim/microscale/BatchedTorchAnalysis.h index 872fe016..55687c7a 100644 --- a/src/mumfim/microscale/BatchedTorchAnalysis.h +++ b/src/mumfim/microscale/BatchedTorchAnalysis.h @@ -6,6 +6,7 @@ #include #include #include +#include #include @@ -25,9 +26,11 @@ namespace mumfim BatchedTorchAnalysis(const std::string & model_path, int nrves) : BatchedRVEAnalysis(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()); 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(); return true; @@ -85,6 +98,10 @@ namespace mumfim KOKKOS_ASSERT(C.extent(0) == this->num_rves_); Kokkos::deep_copy(C.template view(), stiffness_trial_); } + void updateDamageFactor(Kokkos::View 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 F_current_; + Kokkos::View damage_factor_; torch::jit::script::Module model_; diff --git a/src/mumfim/microscale/MultiscaleRVEAnalysis.cc b/src/mumfim/microscale/MultiscaleRVEAnalysis.cc index d5350251..b80fee7c 100644 --- a/src/mumfim/microscale/MultiscaleRVEAnalysis.cc +++ b/src/mumfim/microscale/MultiscaleRVEAnalysis.cc @@ -25,6 +25,21 @@ #endif namespace mumfim { + + static void calculateUpdatedDamageFactor(double failure_stress, double damage_factor, + const Kokkos::View & stress, const Kokkos::View & accepted_damage, + const Kokkos::View & 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 orientation_tensor_normal( "orientation tensor normal", hdrs.size()); - Kokkos::View trial_damage("damage factor", + Kokkos::View trial_damage("damage factor", hdrs.size()); - Kokkos::View accepted_damage("damage factor", + auto trial_damage_h = Kokkos::create_mirror_view(trial_damage); + Kokkos::View 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(); - 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(); material_stiffness.sync(); + 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); diff --git a/src/mumfim/microscale/batched_analysis/BatchedFiberRVEAnalysisExplicitBase.h b/src/mumfim/microscale/batched_analysis/BatchedFiberRVEAnalysisExplicitBase.h index a4a8efb2..acabe573 100644 --- a/src/mumfim/microscale/batched_analysis/BatchedFiberRVEAnalysisExplicitBase.h +++ b/src/mumfim/microscale/batched_analysis/BatchedFiberRVEAnalysisExplicitBase.h @@ -293,5 +293,18 @@ namespace mumfim } } }; + inline Kokkos::View computeVonMisesStress(Kokkos::View stress) { + Kokkos::View 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