From 66b9738c14e61af8a5c6824265f0800161aeb70d Mon Sep 17 00:00:00 2001
From: pabloseleson <selesonpd@ornl.gov>
Date: Fri, 22 Nov 2024 12:16:22 -0500
Subject: [PATCH] Added critical time step for heat transfer

---
 src/CabanaPD_Input.hpp | 87 +++++++++++++++++++++++++++++++++++++++---
 1 file changed, 81 insertions(+), 6 deletions(-)

diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp
index 3105f2e..1ab1897 100644
--- a/src/CabanaPD_Input.hpp
+++ b/src/CabanaPD_Input.hpp
@@ -62,9 +62,13 @@ class Inputs
             inputs["bulk_modulus"]["value"] = K;
         }
 
-        // Check critical time step
+        // Check critical time steps
         // This must be done after the values above are calculated
-        computeCriticalTimeStep();
+        computeCriticalTimeStepMechanics();
+
+        // TODO: Here we need a condition to check whether heat transfer is
+        // included
+        computeCriticalTimeStepHeatTransfer();
 
         int num_steps = tf / dt;
         inputs["num_steps"]["value"] = num_steps;
@@ -177,7 +181,7 @@ class Inputs
             log_err( std::cout, "Units for dx do not match system units." );
     }
 
-    void computeCriticalTimeStep()
+    void computeCriticalTimeStepMechanics()
     {
         // Reference: Silling & Askari, Computers & Structures 83(17–18) (2005):
         // 1526-1535.
@@ -238,11 +242,82 @@ class Inputs
         if ( dt > dt_crit )
         {
             log( std::cout, "WARNING: timestep (", dt,
-                 ") is larger than estimated stable timestep (", dt_crit,
-                 "), using safety factor of ", safety_factor, ".\n" );
+                 ") is larger than estimated stable timestep for mechanics (",
+                 dt_crit, "), using safety factor of ", safety_factor, ".\n" );
+        }
+        // Store in inputs
+        inputs["critical_timestep_mechanics"]["value"] = dt_crit;
+    }
+
+    void computeCriticalTimeStepHeatTransfer()
+    {
+        // The calculation here uses a similar procedure as for the PMB critical
+        // time step calculation
+
+        // Compute particle volume
+        double dx = inputs["dx"]["value"][0];
+        double dy = inputs["dx"]["value"][1];
+        double dz = inputs["dx"]["value"][2];
+        double v_p = dx * dy * dz;
+
+        // Initialize denominator's summation
+        double sum = 0;
+
+        // Run over the neighborhood of a point in the bulk of a body
+        int m = inputs["m"]["value"];
+        double rho = inputs["density"]["value"];
+        double cp = inputs["specific_heat_capacity"]["value"];
+        double delta = inputs["horizon"]["value"];
+
+        for ( int i = -( m + 1 ); i < m + 2; i++ )
+        {
+            // x-component of bond
+            double xi_1 = i * dx;
+
+            for ( int j = -( m + 1 ); j < m + 2; j++ )
+            {
+                // y-component of bond
+                double xi_2 = j * dy;
+
+                for ( int k = -( m + 1 ); k < m + 2; k++ )
+                {
+                    // z-component of bond
+                    double xi_3 = k * dz;
+
+                    // Bond length squared
+                    double r2 = xi_1 * xi_1 + xi_2 * xi_2 + xi_3 * xi_3;
+
+                    // Check if bond is no longer than delta
+                    if ( r2 < delta * delta + 1e-10 )
+                    {
+                        // Check is bond is not 0
+                        if ( r2 > 0 )
+                        {
+                            double xi = std::sqrt( r2 );
+                            double coeff =
+                                model.microconductivity_function( xi );
+
+                            // Compute denominator
+                            sum += v_p * coeff / r2;
+                        }
+                    }
+                }
+            }
+        }
+
+        double safety_factor = inputs["timestep_safety_factor"]["value"];
+        double dt_crit = safety_factor * rho * cp / sum;
+
+        double dt = inputs["timestep"]["value"];
+        if ( dt > dt_crit )
+        {
+            log( std::cout, "WARNING: timestep (", dt,
+                 ") is larger than estimated stable timestep for heat transfer "
+                 "(",
+                 dt_crit, "), using safety factor of ", safety_factor, ".\n" );
         }
         // Store in inputs
-        inputs["critical_timestep"]["value"] = dt_crit;
+        inputs["critical_timestep_heat_transfer"]["value"] = dt_crit;
     }
 
     // Parse JSON file.