Feature neo hookean

.gitignore

@@ -36,3 +36,5 @@ Developer_Tools/sweep_results.json
 build_matrix_*.json
 build_matrix_*.txt
 clangformat.sh
+Tests/**/PreProcess/__pycache__/*
+sweep_results.json

Developer_Tools/Run_All_Tests.py

@@ -302,7 +302,7 @@ def Run_ParameterSweep_1D_Axial_Bar_Vibration(cfg):
         config["use_sycl"] = use_sycl
         config["use_eb"] = use_eb
         config["use_temp"] = use_temp
-
+        config["density"] = 1.0
         config["output_tag"] = output_tag
 
         # 3. Write updated config.json
@@ -527,6 +527,7 @@ def Run_ParameterSweep_1D_HeatConduction(cfg):
         config["use_sycl"] = use_sycl
         config["use_eb"] = use_eb
         config["use_temp"] = use_temp
+        config["density"] = 1.0
 
         # Auto-tag       
         config["output_tag"] = output_tag
@@ -730,6 +731,7 @@ def Run_ParameterSweep_1D_HeatConduction_HeatFlux(cfg):
         config["use_sycl"] = use_sycl
         config["use_eb"] = use_eb
         config["use_temp"] = use_temp
+        config["density"] = 1.0
 
         # Auto-tag       
         config["output_tag"] = output_tag
@@ -936,6 +938,7 @@ def Run_ParameterSweep_1D_HeatConduction_Convective(cfg):
         config["use_sycl"] = use_sycl
         config["use_eb"] = use_eb
         config["use_temp"] = use_temp
+        config["density"] = 1.0
 
         # Auto-tag       
         config["output_tag"] = output_tag
@@ -1145,6 +1148,7 @@ def Run_ParameterSweep_2D_HeatConduction(cfg):
         config["use_temp"] = use_temp
         # Auto-tag
         config["output_tag"] = output_tag
+        config["density"] = 1.0
 
         # 3. Write updated config.json
         with open(os.path.join(test_dir, "./PreProcess/config.json"), "w") as f:
@@ -1355,6 +1359,7 @@ def Run_ParameterSweep_2D_HeatConduction_Cylinder_Dirichlet(cfg):
         config["use_eb"] = use_eb
         config["use_temp"] = use_temp
         config["output_tag"] = output_tag
+        config["density"] = 1.0
 
         # 3. Write updated config.json
         with open(os.path.join(test_dir, "./PreProcess/config.json"), "w") as f:
@@ -1561,6 +1566,7 @@ def Run_ParameterSweep_Dambreak(cfg):
         config["use_sycl"] = use_sycl
         config["use_eb"] = use_eb
         config["use_temp"] = use_temp
+        config["density"] = 997.0
 
         # Auto-tag       
         config["output_tag"] = output_tag
@@ -1768,6 +1774,7 @@ def Run_ParameterSweep_EDC(cfg):
         config["use_temp"] = use_temp
         # Auto-tag       
         config["output_tag"] = output_tag
+        config["density"] = 997.0
 
         # 3. Write updated config.json
         with open(os.path.join(test_dir, "./PreProcess/config.json"), "w") as f:
@@ -2511,7 +2518,7 @@ def _cmake_bool(val):
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["gnumake","cmake"],
+            "build_system": ["cmake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2539,7 +2546,7 @@ def _cmake_bool(val):
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["gnumake","cmake"],
+            "build_system": ["cmake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2564,7 +2571,7 @@ def _cmake_bool(val):
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["gnumake","cmake"],
+            "build_system": ["gnumake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2584,12 +2591,12 @@ def _cmake_bool(val):
         "parameter_space": {
             "dimension": [2],
             "np_per_cell_x": [2],            
-            "order_scheme": [1],            
+            "order_scheme": [2],            
             "stress_update_scheme": ["MUSL"],
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["gnumake","cmake"],
+            "build_system": ["gnumake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2613,7 +2620,7 @@ def _cmake_bool(val):
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["gnumake","cmake"],
+            "build_system": ["gnumake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2663,7 +2670,7 @@ def _cmake_bool(val):
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["cmake","gnumake"],
+            "build_system": ["cmake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2688,7 +2695,7 @@ def _cmake_bool(val):
             "build_with_hdf": [True],
             "output_format": ["hdf5"],
             "filename_prefix": ["mpm_particles"],            
-            "build_system": ["cmake","gnumake"],
+            "build_system": ["cmake"],
             "use_mpi": [True],
             "use_cuda": [False],
             "use_hip": [False],
@@ -2839,16 +2846,16 @@ def _run_parameter_sweeps():
             #Run_ParameterSweep_1D_Axial_Bar_Vibration(cfg)
         elif test_name == "1D_Heat_Conduction":
             print('Nothing to do')
-            #Run_ParameterSweep_1D_HeatConduction(cfg)
+            Run_ParameterSweep_1D_HeatConduction(cfg)
         elif test_name == "1D_Heat_Conduction_HeatFlux":
             print('Nothing to do')
-            #Run_ParameterSweep_1D_HeatConduction_HeatFlux(cfg)
+            Run_ParameterSweep_1D_HeatConduction_HeatFlux(cfg)
         elif test_name == "1D_Heat_Conduction_Convective":
             print('Nothing to do')
-            #Run_ParameterSweep_1D_HeatConduction_Convective(cfg)
+            Run_ParameterSweep_1D_HeatConduction_Convective(cfg)
         elif test_name == "2D_Heat_Conduction":
             print('Nothing to do')
-            #Run_ParameterSweep_2D_HeatConduction(cfg)
+            Run_ParameterSweep_2D_HeatConduction(cfg)
         elif test_name == "2D_Heat_Conduction_Cylinder_Dirichlet":
             print('Nothing to do')
             Run_ParameterSweep_2D_HeatConduction_Cylinder_Dirichlet(cfg)

Source/constitutive_models.H

@@ -37,7 +37,7 @@ linear_elastic(amrex::Real eps[NCOMP_TENSOR],
     // 2D plane stress/strain case
     sigma[XX] = const1 * ((1 - v) * eps[XX] + v * eps[YY]);
     sigma[YY] = const1 * ((1 - v) * eps[YY] + v * eps[XX]);
-    sigma[ZZ] = 0.0; // out-of-plane stress often set to 0 in plane stress
+    sigma[ZZ] = v * (sigma[XX] + sigma[YY]);
     sigma[XY] = const2 * eps[XY];
     sigma[XZ] = sigma[YZ] = 0.0;
 
@@ -84,7 +84,7 @@ linear_elastic_delta(amrex::Real delta_eps[NCOMP_TENSOR],
     // 2D elasticity (plane stress assumption)
     delta_sigma[XX] = const1 * ((1 - v) * delta_eps[XX] + v * delta_eps[YY]);
     delta_sigma[YY] = const1 * ((1 - v) * delta_eps[YY] + v * delta_eps[XX]);
-    delta_sigma[ZZ] = 0.0;
+    delta_sigma[ZZ] = v * (delta_sigma[XX] + delta_sigma[YY]);
     delta_sigma[XY] = const2 * delta_eps[XY];
     delta_sigma[XZ] = delta_sigma[YZ] = 0.0;
 
@@ -136,7 +136,7 @@ Newtonian_Fluid(amrex::Real epsdot[NCOMP_TENSOR],
         2.0 * dyn_visc * (epsdot[XX] - one_by_two * trace_epsdot) - pressure;
     sigma[YY] =
         2.0 * dyn_visc * (epsdot[YY] - one_by_two * trace_epsdot) - pressure;
-    sigma[ZZ] = 0.0;
+    sigma[ZZ] = 2 * dyn_visc * (epsdot[ZZ] - trace_epsdot / 3.0) - pressure;
     sigma[XY] = 2.0 * dyn_visc * epsdot[XY];
     sigma[XZ] = sigma[YZ] = 0.0;
 
@@ -158,4 +158,95 @@ Newtonian_Fluid(amrex::Real epsdot[NCOMP_TENSOR],
 #endif
 }
 
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
+left_cauchy_green(const amrex::Real F[NCOMP_FULLTENSOR],
+                  amrex::Real b[NCOMP_TENSOR])
+{
+    constexpr int d = AMREX_SPACEDIM; // stride into F, NOT 3
+
+    auto row_dot = [&](int i, int j)
+    {
+        amrex::Real s = amrex::Real(0.0);
+        for (int k = 0; k < d; ++k)
+            s += F[i * d + k] * F[j * d + k];
+        return s;
+    };
+
+    b[XX] = row_dot(0, 0);
+#if (AMREX_SPACEDIM >= 2)
+    b[YY] = row_dot(1, 1);
+    b[XY] = row_dot(0, 1);
+#endif
+#if (AMREX_SPACEDIM == 3)
+    b[ZZ] = row_dot(2, 2);
+    b[XZ] = row_dot(0, 2);
+    b[YZ] = row_dot(1, 2);
+#endif
+}
+
+/**
+ * @brief Computes Cauchy stress for a neo-hooken solid
+ *
+ * Implements compressible neo-hookean elasticity in 1D, 2D, or 3D depending on
+ * AMREX_SPACEDIM. Uses Young’s modulus (E) and Poisson’s ratio (v) to form
+ * the constitutive matrix and multiplies it by the strain tensor eps[].
+ *
+ * @param eps    Strain tensor components
+ * @param sigma  Output stress tensor components
+ * @param E      Young’s modulus
+ * @param v      Poisson’s ratio
+ */
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
+neo_hookean(amrex::Real sigma[NCOMP_TENSOR],
+            const amrex::Real F[NCOMP_FULLTENSOR],
+            amrex::Real E,
+            amrex::Real v)
+{
+    amrex::Real mu = E / (two * (eka + v));                     // shear modulus
+    amrex::Real lambda = E * v / ((eka + v) * (eka - two * v)); // first Lamé
+
+    amrex::Real B[NCOMP_TENSOR];
+    left_cauchy_green(F, B); // B = F Fᵀ
+
+    // Jacobian
+    amrex::Real detF = 0.0;
+#if (AMREX_SPACEDIM == 1)
+    detF = F[0];
+#elif (AMREX_SPACEDIM == 2)
+    detF = F[0] * F[3] - F[1] * F[2];
+#else
+    detF = F[0] * (F[4] * F[8] - F[5] * F[7]) -
+           F[1] * (F[3] * F[8] - F[5] * F[6]) +
+           F[2] * (F[3] * F[7] - F[4] * F[6]);
+#endif
+    detF = amrex::max(detF, amrex::Real(1.e-12)); // NaN guard
+
+    const amrex::Real vol = lambda * std::log(detF) / detF; // λ ln(J)/J
+    const amrex::Real muJ = mu / detF;                      // μ/J
+
+#if (AMREX_SPACEDIM == 1)
+    sigma[XX] = muJ * (B[XX] - 1.0) + vol;
+    sigma[YY] = sigma[ZZ] = vol;
+    sigma[XY] = sigma[XZ] = sigma[YZ] = 0.0;
+
+#elif (AMREX_SPACEDIM == 2)
+    // plane strain: F_zz = 1 implicit, so b_zz = 1 → out-of-plane deviatoric =
+    // 0
+    sigma[XX] = muJ * (B[XX] - 1.0) + vol;
+    sigma[YY] = muJ * (B[YY] - 1.0) + vol;
+    sigma[ZZ] = vol; // NOT zero under plane strain
+    sigma[XY] = muJ * B[XY];
+    sigma[XZ] = sigma[YZ] = 0.0;
+
+#elif (AMREX_SPACEDIM == 3)
+    sigma[XX] = muJ * (B[XX] - 1.0) + vol;
+    sigma[YY] = muJ * (B[YY] - 1.0) + vol;
+    sigma[ZZ] = muJ * (B[ZZ] - 1.0) + vol;
+    sigma[XY] = muJ * B[XY];
+    sigma[XZ] = muJ * B[XZ];
+    sigma[YZ] = muJ * B[YZ];
+#endif
+}
+
 #endif

Source/mpm_init.cpp

@@ -178,7 +178,6 @@ void Initialise_Domain(MPMspecs &specs,
     }
     else if (specs.order_scheme == 2)
     {
-        amrex::Print() << "\n Yes the order is 2";
         ng_cells_nodaldata = 3;
 
         // Set directional order-scheme based on periodicity and grid size
@@ -813,6 +812,19 @@ void MPMParticleContainer::InitParticlesFromHDF5(const std::string &filename,
             p.rdata(realData::Dynamic_viscosity) =
                 extra_data.at("Dynamic_viscosity")[local_i];
         }
+        else if (cm_id[local_i] == 2)
+        {
+            p.rdata(realData::E) = extra_data.at("E")[local_i];
+            p.rdata(realData::nu) = extra_data.at("nu")[local_i];
+            p.rdata(realData::Bulk_modulus) = 0.0;
+            p.rdata(realData::Gama_pressure) = 0.0;
+            p.rdata(realData::Dynamic_viscosity) = 0.0;
+        }
+        else
+        {
+            amrex::Abort("\nUnknown constitutive model encountered in "
+                         "InitParticlesFromHDF5.\n");
+        }
 
 #if USE_TEMP
         p.rdata(realData::temperature) = extra_data.at("T")[local_i];
@@ -1066,8 +1078,19 @@ void MPMParticleContainer::InitParticles(const std::string &filename,
                 safe_read(ifs, p.rdata(realData::Dynamic_viscosity),
                           "Error reading Dynamic_viscosity");
             }
+            else if (p.idata(intData::constitutive_model) == 2)
+            {
+                safe_read(ifs, p.rdata(realData::E), "Error reading E");
+                safe_read(ifs, p.rdata(realData::nu), "Error reading nu");
+                p.rdata(realData::Bulk_modulus) = 0.0;
+                p.rdata(realData::Gama_pressure) = 0.0;
+                p.rdata(realData::Dynamic_viscosity) = 0.0;
+            }
             else
             {
+                amrex::Print() << "Error: Constitutive model ID "
+                               << p.idata(intData::constitutive_model)
+                               << " is not recognized.\n";
                 amrex::Abort("Incorrect constitutive model");
             }
 

Source/mpm_particle_container.cpp

@@ -69,6 +69,7 @@ void MPMParticleContainer::apply_constitutive_model(
                     amrex::Real strainrate[NCOMP_TENSOR];
                     amrex::Real strain[NCOMP_TENSOR];
                     amrex::Real stress[NCOMP_TENSOR];
+                    amrex::Real deformation_gradient[NCOMP_FULLTENSOR];
 
                     // Update strain from strainrate
                     for (int d = 0; d < NCOMP_TENSOR; ++d)
@@ -94,6 +95,13 @@ void MPMParticleContainer::apply_constitutive_model(
                         strain[d] = p.rdata(realData::strain + d);
                     }
 
+                    // Get deformation gradient from particle data
+                    for (int comp = 0; comp < NCOMP_FULLTENSOR; ++comp)
+                    {
+                        deformation_gradient[comp] =
+                            p.rdata(realData::deformation_gradient + comp);
+                    }
+
                     if (p.idata(intData::constitutive_model) == 0)
                     {
                         // Elastic solid
@@ -115,6 +123,13 @@ void MPMParticleContainer::apply_constitutive_model(
                                         p.rdata(realData::Dynamic_viscosity),
                                         p.rdata(realData::pressure));
                     }
+                    else if (p.idata(intData::constitutive_model) == 2)
+                    {
+                        // Neo-Hookean solid
+                        neo_hookean(stress, deformation_gradient,
+                                    p.rdata(realData::E),
+                                    p.rdata(realData::nu));
+                    }
 
                     // Write back stress
                     for (int d = 0; d < NCOMP_TENSOR; ++d)
@@ -187,7 +202,7 @@ void MPMParticleContainer::apply_constitutive_model_delta(
 
                 if (p.idata(intData::phase) == 0)
                 {
-                    amrex::Real delta_strain[NCOMP_TENSOR];
+                    amrex::Real delta_strain[NCOMP_TENSOR] = {};
                     amrex::Real delta_stress[NCOMP_TENSOR];
 
                     // Accumulate strain from current strainrate
@@ -240,6 +255,11 @@ void MPMParticleContainer::apply_constitutive_model_delta(
                             "\nDelta strain model for weakly compressible "
                             "fluids not implemented yet.");
                     }
+                    else if (p.idata(intData::constitutive_model) == 2)
+                    {
+                        amrex::Abort("\nDelta strain model for neo hookean "
+                                     "model not implemented yet.");
+                    }
 
                     // Accumulate stress with delta contribution
                     for (int c = 0; c < NCOMP_TENSOR; ++c)