diff --git a/.gitignore b/.gitignore index a2c1b17..d0773ef 100644 --- a/.gitignore +++ b/.gitignore @@ -36,3 +36,5 @@ Developer_Tools/sweep_results.json build_matrix_*.json build_matrix_*.txt clangformat.sh +Tests/**/PreProcess/__pycache__/* +sweep_results.json diff --git a/Developer_Tools/Run_All_Tests.py b/Developer_Tools/Run_All_Tests.py index 635ae77..bbdd3ef 100644 --- a/Developer_Tools/Run_All_Tests.py +++ b/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) diff --git a/Source/constitutive_models.H b/Source/constitutive_models.H index e96b00f..8d58a61 100644 --- a/Source/constitutive_models.H +++ b/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 diff --git a/Source/mpm_init.cpp b/Source/mpm_init.cpp index f375abe..b87be18 100644 --- a/Source/mpm_init.cpp +++ b/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"); } diff --git a/Source/mpm_particle_container.cpp b/Source/mpm_particle_container.cpp index ab6e4ee..9b9ecee 100644 --- a/Source/mpm_particle_container.cpp +++ b/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) diff --git a/Source/mpm_particle_timestep.cpp b/Source/mpm_particle_timestep.cpp index 71a7cd7..c95c72d 100644 --- a/Source/mpm_particle_timestep.cpp +++ b/Source/mpm_particle_timestep.cpp @@ -60,7 +60,8 @@ amrex::Real MPMParticleContainer::Calculate_time_step(MPMspecs &specs) Cs = std::sqrt(p.rdata(realData::Bulk_modulus) / p.rdata(realData::density)); } - else if (p.idata(intData::constitutive_model) == 0) + else if (p.idata(intData::constitutive_model) == 0 or + p.idata(intData::constitutive_model) == 2) { amrex::Real lambda = p.rdata(realData::E) * @@ -72,6 +73,11 @@ amrex::Real MPMParticleContainer::Calculate_time_step(MPMspecs &specs) Cs = std::sqrt((lambda + 2.0 * mu) / p.rdata(realData::density)); } + else + { + amrex::Abort("\nInvalid constitutive model. dt approaching " + "infinity.\n"); + } // Dimension‑aware velocity magnitude amrex::Real velmag = 0.0; diff --git a/Source/mpm_specs.H b/Source/mpm_specs.H index 2b521f8..8a44a46 100644 --- a/Source/mpm_specs.H +++ b/Source/mpm_specs.H @@ -489,15 +489,15 @@ class MPMspecs auto bc_temp_str_to_int = [](const std::string &s) -> int { if (s == "periodic") - return 0; + return BC_TEMP_PERIODIC; if (s == "dirichlet") - return 1; + return BC_TEMP_ISOTHERMAL; if (s == "adiabatic") - return 2; + return BC_TEMP_ADIABATIC; if (s == "heatflux") - return 3; + return BC_TEMP_USERFLUX; if (s == "convective") - return 4; + return BC_TEMP_CONVECTION; amrex::Abort("\nUnknown temperature BC type: " + s); return -1; }; diff --git a/Source/nodal_data_ops.cpp b/Source/nodal_data_ops.cpp index e0f6f97..a0765a8 100644 --- a/Source/nodal_data_ops.cpp +++ b/Source/nodal_data_ops.cpp @@ -857,8 +857,7 @@ void nodal_bcs_temperature(const amrex::Geometry geom, if (bc_type == BC_TEMP_PERIODIC) { - // Periodic. Do nothing - bc_applied = true; + bc_applied = true; // Periodic. Do nothing } else if (bc_type == BC_TEMP_ISOTHERMAL) { @@ -894,9 +893,12 @@ void nodal_bcs_temperature(const amrex::Geometry geom, IntVect nb = nodeid; nb[d] += sign; amrex::Real hc = is_lo ? h_lo_g[d] : h_hi_g[d]; + amrex::Real mk = arr(nodeid, MASS_CONDUCTIVITY); + amrex::Real m = arr(nodeid, MASS_INDEX); amrex::Real Tinf = is_lo ? Tinf_lo_g[d] : Tinf_hi_g[d]; - amrex::Real Bi = hc * dx_g[d]; + amrex::Real k_node = (m > shunya) ? mk / m : eka; + amrex::Real Bi = hc * dx_g[d] / k_node; arr(nodeid, TEMPERATURE) = (arr(nb, TEMPERATURE) + Bi * Tinf) / (1.0 + Bi); bc_applied = true; @@ -1431,7 +1433,10 @@ void apply_udf_nodal_bcs_temperature(const amrex::Geometry &geom, nb[dir] += sign; amrex::Real hc = val0; amrex::Real Tinf = val1; - amrex::Real Bi = hc * dx_g[dir]; + amrex::Real mk = arr(nodeid, MASS_CONDUCTIVITY); + amrex::Real m = arr(nodeid, MASS_INDEX); + amrex::Real k_node = (m > shunya) ? mk / m : eka; + amrex::Real Bi = hc * dx_g[dir] / k_node; arr(nodeid, TEMPERATURE) = (arr(nb, TEMPERATURE) + Bi * Tinf) / (eka + Bi); } diff --git a/Source/utilities.cpp b/Source/utilities.cpp index 0026f98..f504ddd 100644 --- a/Source/utilities.cpp +++ b/Source/utilities.cpp @@ -225,45 +225,18 @@ void Apply_Nodal_BCs_Temperature(amrex::Geometry &geom, amrex::Real t, bool dirichlet_only) { - if (!dirichlet_only) - { - nodal_bcs_temperature( - geom, nodaldata, specs.bclo_temp.data(), specs.bchi_temp.data(), - specs.bc_temp_T_wall_lo.data(), specs.bc_temp_T_wall_hi.data(), - specs.bc_temp_flux_lo.data(), specs.bc_temp_flux_hi.data(), - specs.bc_temp_h_lo.data(), specs.bc_temp_h_hi.data(), - specs.bc_temp_Tinf_lo.data(), specs.bc_temp_Tinf_hi.data()); - compute_udf_temp_at_nodes(geom, specs, t); - apply_udf_nodal_bcs_temperature(geom, nodaldata, specs); -#if USE_EB - if (mpm_ebtools::using_levelset_geometry) - nodal_levelset_bcs_temperature(nodaldata, geom, - /*dirichlet_only=*/false); -#endif - } - else - { - amrex::Vector bclo_dirichlet(AMREX_SPACEDIM, 0); - amrex::Vector bchi_dirichlet(AMREX_SPACEDIM, 0); - for (int d = 0; d < AMREX_SPACEDIM; ++d) - { - bclo_dirichlet[d] = (specs.bclo_temp[d] == 1) ? 1 : 0; - bchi_dirichlet[d] = (specs.bchi_temp[d] == 1) ? 1 : 0; - } - nodal_bcs_temperature( - geom, nodaldata, bclo_dirichlet.data(), bchi_dirichlet.data(), - specs.bc_temp_T_wall_lo.data(), specs.bc_temp_T_wall_hi.data(), - specs.bc_temp_flux_lo.data(), specs.bc_temp_flux_hi.data(), - specs.bc_temp_h_lo.data(), specs.bc_temp_h_hi.data(), - specs.bc_temp_Tinf_lo.data(), specs.bc_temp_Tinf_hi.data()); - compute_udf_temp_at_nodes(geom, specs, t); - apply_udf_nodal_bcs_temperature(geom, nodaldata, specs); + nodal_bcs_temperature( + geom, nodaldata, specs.bclo_temp.data(), specs.bchi_temp.data(), + specs.bc_temp_T_wall_lo.data(), specs.bc_temp_T_wall_hi.data(), + specs.bc_temp_flux_lo.data(), specs.bc_temp_flux_hi.data(), + specs.bc_temp_h_lo.data(), specs.bc_temp_h_hi.data(), + specs.bc_temp_Tinf_lo.data(), specs.bc_temp_Tinf_hi.data()); + compute_udf_temp_at_nodes(geom, specs, t); + apply_udf_nodal_bcs_temperature(geom, nodaldata, specs); #if USE_EB - if (mpm_ebtools::using_levelset_geometry) - nodal_levelset_bcs_temperature(nodaldata, geom, - /*dirichlet_only=*/true); + if (mpm_ebtools::using_levelset_geometry) + nodal_levelset_bcs_temperature(nodaldata, geom, dirichlet_only); #endif - } } #endif @@ -520,7 +493,6 @@ void Initialise_Diagnostic_Streams(MPMspecs &specs) if (specs.do_calculate_mwa_velmag) { - amrex::Print() << "\n Diag vel mag"; std::string fullfilename = specs.diagnostic_output_folder + "/" + specs.file_mwa_velmag; if (amrex::ParallelDescriptor::IOProcessor()) diff --git a/Tests/1D_Axial_Bar_Vibration/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/1D_Axial_Bar_Vibration/PreProcess/Generate_MPs_Inputfile_Generic.py index 910d99f..575c69d 100644 --- a/Tests/1D_Axial_Bar_Vibration/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/1D_Axial_Bar_Vibration/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 1.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1472,14 +783,6 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: diff --git a/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json b/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json index 0269294..919137c 100644 --- a/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json +++ b/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json @@ -38,6 +38,7 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, @@ -86,4 +87,4 @@ "do_calculate_minmaxpos": 0, "write_diag_output_time": 0.01 } -} +} \ No newline at end of file diff --git a/Tests/1D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/1D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py index 910d99f..e75f523 100644 --- a/Tests/1D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/1D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 1.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1472,14 +783,6 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: diff --git a/Tests/1D_Heat_Conduction/PreProcess/config.json b/Tests/1D_Heat_Conduction/PreProcess/config.json index be04500..47a1149 100644 --- a/Tests/1D_Heat_Conduction/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction/PreProcess/config.json @@ -35,6 +35,7 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, @@ -91,4 +92,4 @@ "do_calculate_minmaxpos": 0, "write_diag_output_time": 0.01 } -} +} \ No newline at end of file diff --git a/Tests/1D_Heat_Conduction_Convective/Generate_MPs_and_InputFiles.sh b/Tests/1D_Heat_Conduction_Convective/Generate_MPs_and_InputFiles.sh index 501d658..b856863 100644 --- a/Tests/1D_Heat_Conduction_Convective/Generate_MPs_and_InputFiles.sh +++ b/Tests/1D_Heat_Conduction_Convective/Generate_MPs_and_InputFiles.sh @@ -1 +1 @@ -python3 ./PreProcess/Generate_MPs_Inputfile_Generic.py --config ./PreProcess/config.json +python ./PreProcess/Generate_MPs_Inputfile_Generic.py --config ./PreProcess/config.json diff --git a/Tests/1D_Heat_Conduction_Convective/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/1D_Heat_Conduction_Convective/PreProcess/Generate_MPs_Inputfile_Generic.py index 910d99f..e75f523 100644 --- a/Tests/1D_Heat_Conduction_Convective/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/1D_Heat_Conduction_Convective/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 1.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1472,14 +783,6 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: diff --git a/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json b/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json index 605b645..c8064a3 100644 --- a/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json @@ -35,17 +35,18 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, - "order_scheme": 1, + "order_scheme": 2, "stress_update_scheme": "MUSL", - "output_tag": "1D_Heat_Conduction_Convective_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemcmake_46d367", + "output_tag": "1D_Heat_Conduction_Convective_npcx2_ord2_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemgnumake_9af11f", "input_filename": "Inputs_1DHeatConduction_Convective.inp", "output_format": "hdf5", "materialpoint_filename": "mpm_particles.h5", "plot_to_check": false, "build_with_hdf": true, - "build_system": "cmake", + "build_system": "gnumake", "use_mpi": true, "use_cuda": false, "use_hip": false, @@ -82,6 +83,18 @@ "h": 2.0, "T_inf": 0.0 } + }, + "ylo": { + "mom": "noslip", + "temp": { + "type": "adiabatic" + } + }, + "yhi": { + "mom": "noslip", + "temp": { + "type": "adiabatic" + } } }, "diagnostics": { diff --git a/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/Generate_MPs_Inputfile_Generic.py index 6c25437..e75f523 100644 --- a/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 1.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1426,7 +737,7 @@ def plot_2d(x, y, grid): ax.set_xlabel("x") ax.set_ylabel("y") ax.set_aspect("equal") - plt.show() + #plt.show() @@ -1467,19 +778,11 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_aspect("equal") - plt.show() + #plt.show() # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: diff --git a/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json b/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json index 91aa5ac..bc84798 100644 --- a/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json @@ -35,17 +35,18 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, "stress_update_scheme": "MUSL", - "output_tag": "1D_Heat_Conduction_HeatFlux_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemcmake_9be2c3", + "output_tag": "1D_Heat_Conduction_HeatFlux_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemgnumake_f8e8bc", "input_filename": "Inputs_1DHeatConduction_HeatFlux.inp", "output_format": "hdf5", "materialpoint_filename": "mpm_particles.h5", "plot_to_check": false, "build_with_hdf": true, - "build_system": "cmake", + "build_system": "gnumake", "use_mpi": true, "use_cuda": false, "use_hip": false, diff --git a/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py index 6c25437..e75f523 100644 --- a/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 1.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1426,7 +737,7 @@ def plot_2d(x, y, grid): ax.set_xlabel("x") ax.set_ylabel("y") ax.set_aspect("equal") - plt.show() + #plt.show() @@ -1467,19 +778,11 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_aspect("equal") - plt.show() + #plt.show() # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: diff --git a/Tests/2D_Heat_Conduction/PreProcess/config.json b/Tests/2D_Heat_Conduction/PreProcess/config.json index a0ae488..3363f1d 100644 --- a/Tests/2D_Heat_Conduction/PreProcess/config.json +++ b/Tests/2D_Heat_Conduction/PreProcess/config.json @@ -35,11 +35,12 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, "stress_update_scheme": "MUSL", - "output_tag": "2D_Heat_Conduction_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemcmake_a09896", + "output_tag": "2D_Heat_Conduction_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemgnumake_180bec", "input_filename": "Inputs_2DHeatConduction.inp", "output_format": "ascii", "materialpoint_filename": "mpm_particles.dat", @@ -105,4 +106,4 @@ "do_calculate_minmaxpos": 0, "write_diag_output_time": 0.01 } -} +} \ No newline at end of file diff --git a/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/Generate_MPs_Inputfile_Generic.py index 0aa2e11..5341d92 100644 --- a/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 1.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 1.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1397,7 +708,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1467,19 +778,11 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_aspect("equal") - plt.show() + #plt.show() # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: @@ -2059,8 +1354,7 @@ def temperature_function(x, y, z): dim, data = read_particles_ascii(matpt_filename) x = data["x"] y = data.get("y") # None in 1D - z = data.get("z") # None in 1D/2D - print(data) + z = data.get("z") # None in 1D/2D diff --git a/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/config.json b/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/config.json index ee8664e..bfe0fe2 100644 --- a/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/config.json +++ b/Tests/2D_Heat_Conduction_Cylinder_Dirichlet/PreProcess/config.json @@ -35,6 +35,7 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, diff --git a/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py index 7e82460..5341d92 100644 --- a/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 997.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 997.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 997.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1397,7 +708,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1426,7 +737,7 @@ def plot_2d(x, y, grid): ax.set_xlabel("x") ax.set_ylabel("y") ax.set_aspect("equal") - plt.show() + #plt.show() @@ -1467,19 +778,11 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_aspect("equal") - plt.show() + #plt.show() # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: @@ -2059,8 +1354,7 @@ def temperature_function(x, y, z): dim, data = read_particles_ascii(matpt_filename) x = data["x"] y = data.get("y") # None in 1D - z = data.get("z") # None in 1D/2D - print(data) + z = data.get("z") # None in 1D/2D diff --git a/Tests/Dam_Break/PreProcess/config.json b/Tests/Dam_Break/PreProcess/config.json index 91d5002..5af81f8 100644 --- a/Tests/Dam_Break/PreProcess/config.json +++ b/Tests/Dam_Break/PreProcess/config.json @@ -42,10 +42,11 @@ } } ], + "density": 997.0, "CFL": 0.1, "order_scheme": 1, "stress_update_scheme": "MUSL", - "output_tag": "Dam_Break__dim2_npcx1_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPFalse_BuildSystemgnumake_55e0c2", + "output_tag": "Dam_Break__dim2_npcx1_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPFalse_BuildSystemcmake_4f5515", "input_filename": "Inputs_DamBreak.inp", "output_format": "hdf5", "materialpoint_filename": "mpm_particles.h5", @@ -88,7 +89,7 @@ "write_diag_output_time": 0.01 }, "build_with_hdf": true, - "build_system": "gnumake", + "build_system": "cmake", "use_mpi": true, "use_cuda": false, "use_hip": false, diff --git a/Tests/Dam_Break_With_Obstacles/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/Dam_Break_With_Obstacles/PreProcess/Generate_MPs_Inputfile_Generic.py index 594ab38..e75f523 100644 --- a/Tests/Dam_Break_With_Obstacles/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/Dam_Break_With_Obstacles/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 997.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 997.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 997.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1397,7 +708,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1426,7 +737,7 @@ def plot_2d(x, y, grid): ax.set_xlabel("x") ax.set_ylabel("y") ax.set_aspect("equal") - plt.show() + #plt.show() @@ -1467,19 +778,11 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_aspect("equal") - plt.show() + #plt.show() # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1662,10 +965,7 @@ def write_inputs_file( # Embedded Boundary — Level Sets if ls_bodies: body_names = " ".join(b["name"] for b in ls_bodies) - ls_entries = [ - ("eb2.body_names", body_names), - ("mpm.blevset_output_folder", f'"{output_tag}/"'), - ] + ls_entries = [("eb2.body_names", body_names)] for b in ls_bodies: name = b["name"] ls_entries.append((f"{name}.geom_type", b["geom_type"])) @@ -1883,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1913,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1981,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: @@ -2062,8 +1354,7 @@ def temperature_function(x, y, z): dim, data = read_particles_ascii(matpt_filename) x = data["x"] y = data.get("y") # None in 1D - z = data.get("z") # None in 1D/2D - print(data) + z = data.get("z") # None in 1D/2D diff --git a/Tests/Elastic_disk_collision/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/Elastic_disk_collision/PreProcess/Generate_MPs_Inputfile_Generic.py index 7e82460..575c69d 100644 --- a/Tests/Elastic_disk_collision/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/Elastic_disk_collision/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -181,248 +181,6 @@ def ppc_offsets(N: int) -> np.ndarray: return (2 * i - 1) / (2 * N) -def generate_particles_and_return( - dimensions: int, - grid: dict, - ppc: Tuple[int, ...], - constitutive_model: dict, - enable_temperature: bool, - shape_cfg: Optional[dict], - velocity_function: Callable[[float, float, float], Tuple[float, float, float]], - temperature_function: Optional[ - Callable[[float, float, float], Tuple[float, float, float, float]] - ], - out_particles: str = "mpm_particles.dat", -) -> Tuple[int, float]: - if dimensions not in [1, 2, 3]: - die("dimensions must be 1, 2, or 3") - if len(ppc) != dimensions: - die("ppc tuple length must match dimensions") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - dx1 = dx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin = -0.5 * dx - zmax = 0.5 * dx - nz = 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - - #shape_obj = None if dimensions == 1 else make_shape(shape_cfg, dimensions) - - if shape_cfg is None: - shape_obj = None - else: - shape_obj = make_shape(shape_cfg, dimensions) - - - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - phase = 0 - dens = 997.0 - - cm_type = constitutive_model["type"] - - if cm_type == "elastic": - cm_extra = { - "E": constitutive_model["E"], - "nu": constitutive_model["nu"], - } - cm_id = 0 - elif cm_type == "fluid": - cm_extra = { - "Bulk_modulus": constitutive_model["Bulk_modulus"], - "Gama_pressure": constitutive_model["Gama_pressure"], - "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], - } - cm_id = 1 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - def column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - # velocities - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - npart = 0 - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - npart += 1 - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is None or shape_obj.contains((px, py)): - npart += 1 - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is None or shape_obj.contains((px, py, pz)): - npart += 1 - - with open(out_particles, "w") as f: - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: {npart}\n") - f.write("# " + " ".join(column_names()) + "\n") - - for i in range(nx): - cx = xmin + i * dx - for ox in offsets[0]: - px = cx + ox * dx - - if dimensions == 1: - vx, vy, vz = velocity_function(px, 0.0, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{cm_id:d}" - - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, 0.0, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for j in range(ny): - cy = ymin + j * dy - for oy in offsets[1]: - py = cy + oy * dy - - if dimensions == 2: - if shape_obj is not None and not shape_obj.contains((px, py)): - continue - - vx, vy, vz = velocity_function(px, py, 0.0) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = temperature_function( - px, py, 0.0 - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - continue - - for k in range(nz): - cz = zmin + k * dz - for oz in offsets[2]: - pz = cz + oz * dz - if shape_obj is not None and not shape_obj.contains( - (px, py, pz) - ): - continue - - vx, vy, vz = velocity_function(px, py, pz) - cols = [ - f"{phase:d}", - f"{px:.6e}", - f"{py:.6e}", - f"{pz:.6e}", - f"{rad:.6e}", - f"{dens:.6e}", - f"{vx:.6e}", - f"{vy:.6e}", - f"{vz:.6e}", - f"{cm_id:d}" - ] - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - if enable_temperature: - T, spheat, thermcond, heatsrc = ( - temperature_function(px, py, pz) - ) - cols += [ - f"{T:.6e}", - f"{spheat:.6e}", - f"{thermcond:.6e}", - f"{heatsrc:.6e}", - ] - f.write(" ".join(cols) + "\n") - - print(f"WROTE: {out_particles} with {npart} particles (cm_type={cm_type}, id={cm_id})") - return npart, dx1 - def generate_particle_chunks( dimensions, grid, @@ -430,6 +188,7 @@ def generate_particle_chunks( constitutive_model, enable_temperature, shape_cfg, + density, velocity_function, temperature_function, cell_block=(32, 32, 8), @@ -480,7 +239,7 @@ def generate_particle_chunks( vol_particle = vol_cell / np.prod(ppc) rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 997.0 + dens = density phase = 0 # ------------------------------------------------------------ @@ -501,7 +260,13 @@ def generate_particle_chunks( "Gamma_pressure": constitutive_model["Gamma_pressure"], "Dynamic_viscosity": constitutive_model["Dynamic_viscosity"], } - cm_id = 1 + cm_id = 1 + elif cm_type == "neohookean": + cm_extra = { + "E": constitutive_model["E"], + "nu": constitutive_model["nu"], + } + cm_id = 2 else: # Generic fallback for custom models cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} @@ -783,463 +548,9 @@ def flush_chunk(): if chunk is not None: yield chunk - - -def generate_particles_vectorized( - dimensions, - grid, - ppc, - constitutive_model, - enable_temperature, - shape_cfg, - velocity_function, - temperature_function, - out_particles=None, - output_format="ascii", # "ascii", "hdf5", "memory" - cell_block=(32, 32, 8), -): - import numpy as np - - # Optional import for HDF5 mode - if output_format == "hdf5": - import h5py - - # ------------------------------------------------------------ - # Validate mode - # ------------------------------------------------------------ - if output_format not in ("ascii", "hdf5", "memory"): - raise ValueError("output_format must be 'ascii', 'hdf5', or 'memory'") - - is_ascii = (output_format == "ascii") - is_hdf5 = (output_format == "hdf5") - is_memory = (output_format == "memory") - - if (is_ascii or is_hdf5) and out_particles is None: - raise ValueError("out_particles must be provided for ascii or hdf5 output") - - # ------------------------------------------------------------ - # Grid setup - # ------------------------------------------------------------ - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - if dimensions >= 2: - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - dy = (ymax - ymin) / ny - else: - ymin, ymax, ny = 0.0, 1.0, 1 - dy = 1.0 - - if dimensions == 3: - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - dz = (zmax - zmin) / nz - else: - zmin, zmax, nz = -0.5 * dx, 0.5 * dx, 1 - dz = (zmax - zmin) / nz - - offsets = [ppc_offsets(ppc[d]) for d in range(dimensions)] - shape_obj = None if shape_cfg is None else make_shape(shape_cfg, dimensions) - - # ------------------------------------------------------------ - # Volume, radius, density - # ------------------------------------------------------------ - if dimensions == 1: - vol_cell = dx - vol_particle = dx / ppc[0] - elif dimensions == 2: - vol_cell = dx * dy - vol_particle = vol_cell / (ppc[0] * ppc[1]) - else: - vol_cell = dx * dy * dz - vol_particle = vol_cell / np.prod(ppc) - - rad = (3.0 / 4.0 * vol_particle / np.pi) ** (1.0 / 3.0) - dens = 997.0 - phase = 0 - - # ------------------------------------------------------------ - # Constitutive model - # ------------------------------------------------------------ - cm_type = constitutive_model["type"] - if cm_type == "elastic": - cm_extra = {"E": constitutive_model["E"], "nu": constitutive_model["nu"]} - cm_id = 0 - else: - cm_extra = {k: v for k, v in constitutive_model.items() if k != "type"} - cm_id = -1 - - # ------------------------------------------------------------ - # ASCII header helper - # ------------------------------------------------------------ - def ascii_column_names(): - cols = ["phase", "x"] - if dimensions >= 2: - cols.append("y") - if dimensions == 3: - cols.append("z") - - if dimensions == 1: - cols += ["vx"] - elif dimensions == 2: - cols += ["vx", "vy"] - else: - cols += ["vx", "vy", "vz"] - - cols += ["radius", "density"] - cols += list(cm_extra.keys()) - - if enable_temperature: - cols += ["T", "spheat", "thermcond", "heatsrc"] - - return cols - - # ------------------------------------------------------------ - # Output mode setup - # ------------------------------------------------------------ - if is_memory: - mem = {k: [] for k in ["x","y","z","vx","vy","vz","radius","density","cm_id"]} - mem["phase"] = [] - for k in cm_extra.keys(): - mem[k] = [] - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - mem[k] = [] - - elif is_ascii: - f = open(out_particles, "w") - # We do NOT know npart yet → write placeholder, fix later - f.write(f"dim: {dimensions}\n") - f.write(f"number_of_material_points: 0\n") - f.write("# " + " ".join(ascii_column_names()) + "\n") - ascii_count = 0 - - elif is_hdf5: - h5 = h5py.File(out_particles, "w") - h5["dim"] = dimensions - h5["number_of_material_points"] = 0 - - def create_dset(name): - return h5.create_dataset(name, shape=(0,), maxshape=(None,), dtype="f8") - - dsets = {} - for name in ["phase","x","radius","density","vx","cm_id"]: - dsets[name] = create_dset(name) - if dimensions >= 2: - dsets["y"] = create_dset("y") - dsets["vy"] = create_dset("vy") - if dimensions == 3: - dsets["z"] = create_dset("z") - dsets["vz"] = create_dset("vz") - - for k in cm_extra.keys(): - dsets[k] = create_dset(k) - - if enable_temperature: - for k in ["T","spheat","thermcond","heatsrc"]: - dsets[k] = create_dset(k) - - buf = {k: [] for k in dsets.keys()} - total_npart = 0 - - def flush(): - nonlocal total_npart - n = len(buf["x"]) - if n == 0: - return - old = total_npart - new = old + n - for name, dset in dsets.items(): - dset.resize((new,)) - dset[old:new] = np.asarray(buf[name]) - buf[name].clear() - total_npart = new - - # ------------------------------------------------------------ - # Vectorized block generator (2D only) - # ------------------------------------------------------------ - def block_2d(ix0, ix1, iy0, iy1): - ix = np.arange(ix0, ix1) - iy = np.arange(iy0, iy1) - cx = xmin + ix * dx - cy = ymin + iy * dy - CX, CY = np.meshgrid(cx, cy, indexing="ij") - PX = CX[:, :, None] + offsets[0][None, None, :] * dx - PY = CY[:, :, None] + offsets[1][None, None, :] * dy - PX = PX.ravel() - PY = PY.ravel() - if shape_obj is not None: - mask = np.array([shape_obj.contains((x, y)) for x, y in zip(PX, PY)]) - PX = PX[mask] - PY = PY[mask] - return PX, PY - - # ------------------------------------------------------------ - # Main loop - # ------------------------------------------------------------ - bx, by, bz = cell_block - - for ix0 in range(0, nx, bx): - ix1 = min(ix0 + bx, nx) - - for iy0 in range(0, ny, by): - iy1 = min(iy0 + by, ny) - - PX, PY = block_2d(ix0, ix1, iy0, iy1) - PZ = np.zeros_like(PX) - - for px, py, pz in zip(PX, PY, PZ): - vx, vy, vz = velocity_function(px, py, pz) - - if enable_temperature: - T0, SP0, K0, Q0 = temperature_function(px, py, pz) - - # ------------------------- - # MEMORY MODE - # ------------------------- - if is_memory: - mem["phase"].append(phase) - mem["x"].append(px) - mem["y"].append(py) - mem["z"].append(pz) - mem["vx"].append(vx) - mem["vy"].append(vy) - mem["vz"].append(vz) - mem["radius"].append(rad) - mem["density"].append(dens) - mem["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - mem[k].append(v) - if enable_temperature: - mem["T"].append(T0) - mem["spheat"].append(SP0) - mem["thermcond"].append(K0) - mem["heatsrc"].append(Q0) - - # ------------------------- - # ASCII MODE - # ------------------------- - elif is_ascii: - cols = [f"{phase:d}", f"{px:.6e}"] - if dimensions >= 2: - cols.append(f"{py:.6e}") - if dimensions == 3: - cols.append(f"{pz:.6e}") - - if dimensions == 1: - cols.append(f"{vx:.6e}") - elif dimensions == 2: - cols += [f"{vx:.6e}", f"{vy:.6e}"] - else: - cols += [f"{vx:.6e}", f"{vy:.6e}", f"{vz:.6e}"] - - cols += [f"{rad:.6e}", f"{dens:.6e}"] - - for v in cm_extra.values(): - cols.append(f"{v:.6e}") - - if enable_temperature: - cols += [ - f"{T0:.6e}", - f"{SP0:.6e}", - f"{K0:.6e}", - f"{Q0:.6e}", - ] - - f.write(" ".join(cols) + "\n") - ascii_count += 1 - - # ------------------------- - # HDF5 MODE - # ------------------------- - elif is_hdf5: - buf["phase"].append(phase) - buf["x"].append(px) - if dimensions >= 2: - buf["y"].append(py) - if dimensions == 3: - buf["z"].append(pz) - buf["vx"].append(vx) - if dimensions >= 2: - buf["vy"].append(vy) - if dimensions == 3: - buf["vz"].append(vz) - buf["radius"].append(rad) - buf["density"].append(dens) - buf["cm_id"].append(cm_id) - for k, v in cm_extra.items(): - buf[k].append(v) - if enable_temperature: - buf["T"].append(T0) - buf["spheat"].append(SP0) - buf["thermcond"].append(K0) - buf["heatsrc"].append(Q0) - - if is_hdf5: - flush() - - # ------------------------------------------------------------ - # Finalize - # ------------------------------------------------------------ - if is_memory: - return {k: np.asarray(v) for k, v in mem.items()} - - elif is_ascii: - f.close() - # Fix header count - with open(out_particles, "r+") as f2: - lines = f2.readlines() - lines[1] = f"number_of_material_points: {ascii_count}\n" - f2.seek(0) - f2.writelines(lines) - return ascii_count, dx - - elif is_hdf5: - h5["number_of_material_points"][...] = total_npart - h5.close() - return total_npart - # ------------------------------------------------------------ # Plotting # ------------------------------------------------------------ -def plot_material_points( - points: np.ndarray, - grid: dict, - dimensions: int, - output_tag: str, - *, - slice_axis: Optional[str] = None, - slice_value: Optional[float] = None, - figsize=(8, 6), -): - fig, ax = plt.subplots(figsize=figsize) - - if dimensions == 1: - x = points[:, 0] - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - dx = (xmax - xmin) / nx - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - - ax.axvline(xmin, color="black", linewidth=2.5) - ax.axvline(xmax, color="black", linewidth=2.5) - - ax.plot(x, np.zeros_like(x), "o", markersize=4) - ax.set_ylim(-0.1, 0.1) - ax.set_xlabel("x") - ax.set_title("1D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 2: - x = points[:, 0] - y = points[:, 1] - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - - for i in range(nx + 1): - ax.axvline(xmin + i * dx, color="lightgray", linewidth=0.8) - for j in range(ny + 1): - ax.axhline(ymin + j * dy, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin, ymin), - xmax - xmin, - ymax - ymin, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x, y, "o", markersize=3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_aspect("equal") - ax.set_title("2D Material Points with Grid + Boundary") - plt.savefig(output_tag) - return - - if dimensions == 3: - if slice_axis not in ["x", "y", "z"]: - die("For 3D visualization, slice_axis must be 'x', 'y', or 'z'") - if slice_value is None: - die("For 3D visualization, slice_value must be provided") - - xmin, xmax, nx = grid["xmin"], grid["xmax"], grid["nx"] - ymin, ymax, ny = grid["ymin"], grid["ymax"], grid["ny"] - zmin, zmax, nz = grid["zmin"], grid["zmax"], grid["nz"] - - dx = (xmax - xmin) / nx - dy = (ymax - ymin) / ny - dz = (zmax - zmin) / nz - - if slice_axis == "x": - lo = slice_value - dx - hi = slice_value + dx - mask = (points[:, 0] >= lo) & (points[:, 0] <= hi) - pts = points[mask] - x2 = pts[:, 1] - y2 = pts[:, 2] - xlabel, ylabel = "y", "z" - xmin2, xmax2 = ymin, ymax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dy, dz - nx2, ny2 = ny, nz - - elif slice_axis == "y": - lo = slice_value - dy - hi = slice_value + dy - mask = (points[:, 1] >= lo) & (points[:, 1] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 2] - xlabel, ylabel = "x", "z" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = zmin, zmax - dx2, dy2 = dx, dz - nx2, ny2 = nx, nz - - else: - lo = slice_value - dz - hi = slice_value + dz - mask = (points[:, 2] >= lo) & (points[:, 2] <= hi) - pts = points[mask] - x2 = pts[:, 0] - y2 = pts[:, 1] - xlabel, ylabel = "x", "y" - xmin2, xmax2 = xmin, xmax - ymin2, ymax2 = ymin, ymax - dx2, dy2 = dx, dy - nx2, ny2 = nx, ny - - for i in range(nx2 + 1): - ax.axvline(xmin2 + i * dx2, color="lightgray", linewidth=0.8) - for j in range(ny2 + 1): - ax.axhline(ymin2 + j * dy2, color="lightgray", linewidth=0.8) - - rect = patches.Rectangle( - (xmin2, ymin2), - xmax2 - xmin2, - ymax2 - ymin2, - linewidth=2.5, - edgecolor="black", - facecolor="none", - ) - ax.add_patch(rect) - - ax.plot(x2, y2, "o", markersize=3) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_aspect("equal") - ax.set_title(f"3D Slice at {slice_axis}={slice_value} (±1 cell) with Boundary") - plt.savefig(output_tag) - return - def read_grid_from_input(filename): grid = {} @@ -1397,7 +708,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1426,7 +737,7 @@ def plot_2d(x, y, grid): ax.set_xlabel("x") ax.set_ylabel("y") ax.set_aspect("equal") - plt.show() + #plt.show() @@ -1467,19 +778,11 @@ def plot_3d_slice(x, y, z, grid, slice_axis="z", slice_value=None): ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_aspect("equal") - plt.show() + #plt.show() # ------------------------------------------------------------ # Helpers: load particles, write inputs, auto-tag # ------------------------------------------------------------ -def load_particle_positions(filename: str, dimensions: int) -> np.ndarray: - if dimensions == 1: - pts = np.loadtxt(filename, comments="#", skiprows=3, usecols=[1]) - return pts.reshape(-1, 1) - if dimensions == 2: - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2]) - return np.loadtxt(filename, comments="#", skiprows=3, usecols=[1, 2, 3]) - def write_block(f, entries, comment=None): """ @@ -1880,18 +1183,6 @@ def write_particles_hdf5_streaming(filename, chunk_iter, dimensions): return total -def write_particles_hdf5(filename, particles): - import h5py - import numpy as np - - with h5py.File(filename, "w") as h5: - h5["dim"] = 2 - h5["number_of_material_points"] = len(particles["x"]) - - for key, arr in particles.items(): - h5.create_dataset(key, data=np.asarray(arr)) - - # ------------------------------------------------------------ # Main # ------------------------------------------------------------ @@ -1910,10 +1201,13 @@ def main(): order_scheme = cfg["order_scheme"] stress_update_scheme = cfg["stress_update_scheme"] output_tag = cfg.get("output_tag", "").strip() + if(output_tag==""): + output_tag = make_auto_tag_from_cfg(cfg) input_filename = cfg["input_filename"] matpt_filename = cfg["materialpoint_filename"] plot_to_check = cfg["plot_to_check"] CFL = cfg["CFL"] + density = cfg["density"] # user choice: "ascii" or "hdf5" output_format = cfg.get("output_format", "hdf5").lower() @@ -1978,6 +1272,7 @@ def temperature_function(x, y, z): constitutive_model=cm_cfg, enable_temperature=enable_temperature, shape_cfg=shape_cfg, + density=density, velocity_function=velocity_function, temperature_function=temperature_function, # you can tune these if needed: @@ -2059,8 +1354,7 @@ def temperature_function(x, y, z): dim, data = read_particles_ascii(matpt_filename) x = data["x"] y = data.get("y") # None in 1D - z = data.get("z") # None in 1D/2D - print(data) + z = data.get("z") # None in 1D/2D diff --git a/Tests/Elastic_disk_collision/PreProcess/config.json b/Tests/Elastic_disk_collision/PreProcess/config.json index bbd8958..b04558e 100644 --- a/Tests/Elastic_disk_collision/PreProcess/config.json +++ b/Tests/Elastic_disk_collision/PreProcess/config.json @@ -70,10 +70,11 @@ } } ], + "density": 997.0, "CFL": 0.1, "order_scheme": 1, "stress_update_scheme": "MUSL", - "output_tag": "Elastic_disk_collision__dim2_npcx4_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPFalse_BuildSystemgnumake_34471f", + "output_tag": "Elastic_disk_collision__dim2_npcx4_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPFalse_BuildSystemcmake_22c3e0", "input_filename": "Inputs_ElasticDiskCollision.inp", "output_format": "hdf5", "materialpoint_filename": "mpm_particles.h5", @@ -93,6 +94,7 @@ 0.0, 0.0 ], + "density": 1.0, "boundary_conditions": { "xlo": { "mom": "periodic" @@ -116,7 +118,7 @@ "write_diag_output_time": 0.01 }, "build_with_hdf": true, - "build_system": "gnumake", + "build_system": "cmake", "use_mpi": true, "use_cuda": false, "use_hip": false,