From 1a2bf8e96d784ceaf7f00cada214c82767f49a00 Mon Sep 17 00:00:00 2001 From: SreejithNREL Date: Tue, 11 Aug 2026 20:48:39 -0600 Subject: [PATCH 01/11] corrected bug about temperature bcs. Earlier periodic and adiabatic bcs where considered in the same code block. a clear bug. corrected it. added constants to make bc type more clear and coder-friendly --- .gitignore | 2 ++ Source/constants.H | 6 ++++++ Source/nodal_data_ops.cpp | 19 ++++++++++++------- 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index c6885e6..a2c1b17 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ .project .pydevproject +.cproject .DS_Store .clang-format @@ -34,3 +35,4 @@ Developer_Tools/sweep_results.json build_matrix_*.json build_matrix_*.txt +clangformat.sh diff --git a/Source/constants.H b/Source/constants.H index 77354bf..db37469 100644 --- a/Source/constants.H +++ b/Source/constants.H @@ -92,6 +92,12 @@ #define BC_PARTIALSLIPWALL 3 #define BC_OUTFLOW 4 +#define BC_TEMP_PERIODIC 0 +#define BC_TEMP_ISOTHERMAL 1 +#define BC_TEMP_ADIABATIC 2 +#define BC_TEMP_USERFLUX 3 +#define BC_TEMP_CONVECTION 4 + #define USL 0 #define MUSL 1 diff --git a/Source/nodal_data_ops.cpp b/Source/nodal_data_ops.cpp index 6b7ed48..e0f6f97 100644 --- a/Source/nodal_data_ops.cpp +++ b/Source/nodal_data_ops.cpp @@ -855,14 +855,19 @@ void nodal_bcs_temperature(const amrex::Geometry geom, int bc_type = is_lo ? bclo[d] : bchi[d]; int sign = is_lo ? 1 : -1; - if (bc_type == 1) + if (bc_type == BC_TEMP_PERIODIC) + { + // Periodic. Do nothing + bc_applied = true; + } + else if (bc_type == BC_TEMP_ISOTHERMAL) { amrex::Real Tw = is_lo ? T_wall_lo_g[d] : T_wall_hi_g[d]; arr(nodeid, TEMPERATURE) = Tw; bc_applied = true; } - else if (bc_type == 2 || bc_type == 0) + else if (bc_type == BC_TEMP_ADIABATIC) { if (!bc_applied) { @@ -872,7 +877,7 @@ void nodal_bcs_temperature(const amrex::Geometry geom, bc_applied = true; } } - else if (bc_type == 3) + else if (bc_type == BC_TEMP_USERFLUX) { IntVect nb = nodeid; nb[d] += sign; @@ -884,7 +889,7 @@ void nodal_bcs_temperature(const amrex::Geometry geom, arr(nb, TEMPERATURE) + q * dx_g[d] / k_node; bc_applied = true; } - else if (bc_type == 4) + else if (bc_type == BC_TEMP_CONVECTION) { IntVect nb = nodeid; nb[d] += sign; @@ -1406,11 +1411,11 @@ void apply_udf_nodal_bcs_temperature(const amrex::Geometry &geom, amrex::Real val0 = udf_ptr[flat]; amrex::Real val1 = udf_ptr[flat + 1]; - if (bc_type == 1) + if (bc_type == BC_TEMP_ISOTHERMAL) { arr(nodeid, TEMPERATURE) = val0; } - else if (bc_type == 3) + else if (bc_type == BC_TEMP_USERFLUX) { IntVect nb = nodeid; nb[dir] += sign; @@ -1420,7 +1425,7 @@ void apply_udf_nodal_bcs_temperature(const amrex::Geometry &geom, arr(nodeid, TEMPERATURE) = arr(nb, TEMPERATURE) + val0 * dx_g[dir] / k_node; } - else if (bc_type == 4) + else if (bc_type == BC_TEMP_CONVECTION) { IntVect nb = nodeid; nb[dir] += sign; From 341a76ec4ea83cfa8a815f92a05bae6417cae46b Mon Sep 17 00:00:00 2001 From: SreejithNREL Date: Tue, 11 Aug 2026 21:07:50 -0600 Subject: [PATCH 02/11] added neohookean constitutive model --- Source/constitutive_models.H | 97 ++++++++++++++++++++++++++++++- Source/mpm_init.cpp | 11 ++++ Source/mpm_particle_container.cpp | 15 +++++ Source/mpm_particle_timestep.cpp | 3 +- 4 files changed, 122 insertions(+), 4 deletions(-) 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..db763f4 100644 --- a/Source/mpm_init.cpp +++ b/Source/mpm_init.cpp @@ -1066,8 +1066,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..99082fb 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) diff --git a/Source/mpm_particle_timestep.cpp b/Source/mpm_particle_timestep.cpp index 71a7cd7..b36b192 100644 --- a/Source/mpm_particle_timestep.cpp +++ b/Source/mpm_particle_timestep.cpp @@ -55,7 +55,8 @@ amrex::Real MPMParticleContainer::Calculate_time_step(MPMspecs &specs) if (p.idata(intData::phase) == 0) { amrex::Real Cs = 0.0; - if (p.idata(intData::constitutive_model) == 1) + if (p.idata(intData::constitutive_model) == 1 or + p.idata(intData::constitutive_model) == 2) { Cs = std::sqrt(p.rdata(realData::Bulk_modulus) / p.rdata(realData::density)); From 0ce697f8feb85d1321171af74c121e8f1d2f23fd Mon Sep 17 00:00:00 2001 From: SreejithNREL Date: Wed, 12 Aug 2026 15:20:51 -0600 Subject: [PATCH 03/11] I am testing this branch before raising a PR. The Test cases 1D_HeatConduction_Convection and Adiabatic are running into bugs. I need to take a look at it. Commiting this dirty version --- Developer_Tools/Run_All_Tests.py | 19 +- Source/mpm_specs.H | 13 +- Source/nodal_data_ops.cpp | 23 +- .../Generate_MPs_Inputfile_Generic.py | 723 +---------------- .../PreProcess/config.json | 9 +- .../Generate_MPs_Inputfile_Generic.py | 723 +---------------- .../1D_Heat_Conduction/PreProcess/config.json | 9 +- .../Generate_MPs_and_InputFiles.sh | 2 +- .../Generate_MPs_Inputfile_Generic.py | 723 +---------------- .../PreProcess/config.json | 23 +- .../Generate_MPs_Inputfile_Generic.py | 727 +---------------- .../PreProcess/config.json | 1 + .../Generate_MPs_Inputfile_Generic.py | 727 +---------------- .../2D_Heat_Conduction/PreProcess/config.json | 9 +- .../Generate_MPs_Inputfile_Generic.py | 730 +---------------- .../PreProcess/config.json | 1 + .../Generate_MPs_Inputfile_Generic.py | 732 +---------------- Tests/Dam_Break/PreProcess/config.json | 1 + .../Generate_MPs_Inputfile_Generic.py | 737 +----------------- .../Generate_MPs_Inputfile_Generic.py | 732 +---------------- .../PreProcess/config.json | 5 +- 21 files changed, 154 insertions(+), 6515 deletions(-) diff --git a/Developer_Tools/Run_All_Tests.py b/Developer_Tools/Run_All_Tests.py index 635ae77..df6def7 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: @@ -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], @@ -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], @@ -2851,13 +2858,13 @@ def _run_parameter_sweeps(): #Run_ParameterSweep_2D_HeatConduction(cfg) elif test_name == "2D_Heat_Conduction_Cylinder_Dirichlet": print('Nothing to do') - Run_ParameterSweep_2D_HeatConduction_Cylinder_Dirichlet(cfg) + #Run_ParameterSweep_2D_HeatConduction_Cylinder_Dirichlet(cfg) elif test_name == "Dam_Break": print('Nothing to do') #Run_ParameterSweep_Dambreak(cfg) elif test_name == "Elastic_disk_collision": print('Nothing to do') - #Run_ParameterSweep_EDC(cfg) + Run_ParameterSweep_EDC(cfg) # Save results _sweep_results_path = os.path.join(ROOT, "sweep_results.json") diff --git a/Source/mpm_specs.H b/Source/mpm_specs.H index 2b521f8..2d13fae 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; }; @@ -517,10 +517,13 @@ class MPMspecs if (pp.query(face_lo_temp_keys[d].c_str(), bc_temp_str)) { bclo_temp[d] = bc_temp_str_to_int(bc_temp_str); + amrex::Print()<<"\n bc low string = "< 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +777,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 +1177,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 +1195,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 +1266,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..e2261f8 100644 --- a/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json +++ b/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json @@ -38,15 +38,16 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, "stress_update_scheme": "MUSL", "output_tag": "1D_Axial_Bar_Vibration_dim1_npcx2_ord1_flip1.0_susMUSL_CFL0.1_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPFalse_BuildSystemcmake_44edc4", "input_filename": "Inputs_1DAxialBarVibration.inp", - "build_with_hdf": false, - "output_format": "ascii", - "materialpoint_filename": "mpm_particles.dat", + "build_with_hdf": true, + "output_format": "hdf5", + "materialpoint_filename": "mpm_particles.h5", "plot_to_check": true, "build_system": "cmake", "use_mpi": true, @@ -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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +777,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 +1177,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 +1195,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 +1266,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..cde3982 100644 --- a/Tests/1D_Heat_Conduction/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction/PreProcess/config.json @@ -35,16 +35,17 @@ } } ], + "density": 1.0, "CFL": 0.1, "alpha_pic_flip": 1.0, "order_scheme": 1, "stress_update_scheme": "MUSL", "output_tag": "1D_Heat_Conduction_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemcmake_868edf", "input_filename": "Inputs_1DHeatConduction.inp", - "output_format": "ascii", - "materialpoint_filename": "mpm_particles.dat", + "output_format": "hdf5", + "materialpoint_filename": "mpm_particles.h5", "plot_to_check": true, - "build_with_hdf": false, + "build_with_hdf": true, "build_system": "cmake", "use_mpi": true, "use_cuda": false, @@ -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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +777,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 +1177,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 +1195,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 +1266,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..7fdca15 100644 --- a/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json @@ -9,8 +9,8 @@ "ny": 4 }, "ppc": [ - 2, - 2 + 1, + 1 ], "bodies": [ { @@ -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_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": { @@ -92,4 +105,4 @@ "do_calculate_minmaxpos": 0, "write_diag_output_time": 0.01 } -} \ No newline at end of file +} 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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +731,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 +772,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 +1177,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 +1195,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 +1266,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..090696e 100644 --- a/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction_HeatFlux/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/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py index 6c25437..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +731,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 +772,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 +1177,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 +1195,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 +1266,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..775a880 100644 --- a/Tests/2D_Heat_Conduction/PreProcess/config.json +++ b/Tests/2D_Heat_Conduction/PreProcess/config.json @@ -35,16 +35,17 @@ } } ], + "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", "input_filename": "Inputs_2DHeatConduction.inp", - "output_format": "ascii", - "materialpoint_filename": "mpm_particles.dat", + "output_format": "hdf5", + "materialpoint_filename": "mpm_particles.h5", "plot_to_check": false, - "build_with_hdf": false, + "build_with_hdf": true, "build_system": "cmake", "use_mpi": true, "use_cuda": false, @@ -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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +702,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1467,19 +772,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 +1177,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 +1195,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 +1266,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 +1348,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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +702,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1426,7 +731,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 +772,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 +1177,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 +1195,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 +1266,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 +1348,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..a50a190 100644 --- a/Tests/Dam_Break/PreProcess/config.json +++ b/Tests/Dam_Break/PreProcess/config.json @@ -42,6 +42,7 @@ } } ], + "density": 997.0, "CFL": 0.1, "order_scheme": 1, "stress_update_scheme": "MUSL", 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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +702,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1426,7 +731,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 +772,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 +959,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 +1177,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 +1195,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 +1266,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 +1348,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..e480a38 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 # ------------------------------------------------------------ @@ -783,463 +542,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 +702,7 @@ def plot_1d(x, grid): ax.grid(False) plt.tight_layout() - plt.show() + #plt.show() # ------------------------------------------------------------ @@ -1426,7 +731,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 +772,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 +1177,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 +1195,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 +1266,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 +1348,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..521196e 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", @@ -116,7 +117,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, From 7c9c1d0b7c4a2c9ecacbf56d209e1300deb5e4d4 Mon Sep 17 00:00:00 2001 From: "Sreejith N A (NLR)" <98907926+SreejithNREL@users.noreply.github.com> Date: Tue, 11 Aug 2026 21:24:39 -0600 Subject: [PATCH 04/11] corrected bug about temperature bcs. Earlier periodic and adiabatic bcs where considered in the same code block. a clear bug. corrected it. added constants to make bc type more clear and coder-friendly (#19) Co-authored-by: SreejithNREL --- Source/nodal_data_ops.cpp | 23 ++++------------------- 1 file changed, 4 insertions(+), 19 deletions(-) diff --git a/Source/nodal_data_ops.cpp b/Source/nodal_data_ops.cpp index 5e83126..e0f6f97 100644 --- a/Source/nodal_data_ops.cpp +++ b/Source/nodal_data_ops.cpp @@ -855,20 +855,13 @@ void nodal_bcs_temperature(const amrex::Geometry geom, int bc_type = is_lo ? bclo[d] : bchi[d]; int sign = is_lo ? 1 : -1; - if(d==0 and is_hi==true) - { - //amrex::Print()<<"\n At d = "< Date: Tue, 11 Aug 2026 22:01:26 -0600 Subject: [PATCH 05/11] Bug fix splines periodic (#20) * corrected bug about temperature bcs. Earlier periodic and adiabatic bcs where considered in the same code block. a clear bug. corrected it. added constants to make bc type more clear and coder-friendly * corrected a bug where in splines were using a wrong stencil for periodic bcs --------- Co-authored-by: SreejithNREL --- Source/interpolants.H | 59 +++++++++++++++++++------------- Source/mpm_particle_grid_ops.cpp | 28 ++++++++------- 2 files changed, 51 insertions(+), 36 deletions(-) diff --git a/Source/interpolants.H b/Source/interpolants.H index 384c7d5..f02507b 100644 --- a/Source/interpolants.H +++ b/Source/interpolants.H @@ -685,6 +685,7 @@ quadratic_interp(amrex::Real xp[AMREX_SPACEDIM], const amrex::GpuArray dx, amrex::Array4 nodaldata, int comp, + const amrex::GpuArray &periodic, const int *lo, const int *hi) { @@ -712,7 +713,7 @@ quadratic_interp(amrex::Real xp[AMREX_SPACEDIM], else shapetypex = 3; rx = (xp[XDIR] - (plo[XDIR] + (i + l) * dx[XDIR])) / dx[XDIR]; - lval = quadraticspline_1d(shapetypex, rx); + lval = quadraticspline_1d(periodic[0] ? 3 : shapetypex, rx); value += lval * nodaldata(i + l, 0, 0, comp); } @@ -757,8 +758,8 @@ quadratic_interp(amrex::Real xp[AMREX_SPACEDIM], rx = (xp[XDIR] - (plo[XDIR] + (i + l) * dx[XDIR])) / dx[XDIR]; ry = (xp[YDIR] - (plo[YDIR] + (j + m) * dx[YDIR])) / dx[YDIR]; - lval = quadraticspline_1d(shapetypex, rx); - mval = quadraticspline_1d(shapetypey, ry); + lval = quadraticspline_1d(periodic[0] ? 3 : shapetypex, rx); + mval = quadraticspline_1d(periodic[1] ? 3 : shapetypey, ry); value += lval * mval * nodaldata(i + l, j + m, 0, comp); } @@ -824,9 +825,9 @@ quadratic_interp(amrex::Real xp[AMREX_SPACEDIM], ry = (xp[YDIR] - (plo[YDIR] + (j + m) * dx[YDIR])) / dx[YDIR]; rz = (xp[ZDIR] - (plo[ZDIR] + (k + n) * dx[ZDIR])) / dx[ZDIR]; - lval = quadraticspline_1d(shapetypex, rx); - mval = quadraticspline_1d(shapetypey, ry); - nval = quadraticspline_1d(shapetypez, rz); + lval = quadraticspline_1d(periodic[0] ? 3 : shapetypex, rx); + mval = quadraticspline_1d(periodic[1] ? 3 : shapetypey, ry); + nval = quadraticspline_1d(periodic[2] ? 3 : shapetypez, rz); value += lval * mval * nval * nodaldata(i + l, j + m, k + n, comp); } @@ -869,6 +870,7 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], amrex::Array4 nodaldata, int comp0, int comp1, + const amrex::GpuArray &periodic, const int *lo, const int *hi, amrex::Real vals[2]) @@ -893,7 +895,7 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], shapetypex = 3; amrex::Real rx = (xp[XDIR] - (plo[XDIR] + (i + l) * dx[XDIR])) / dx[XDIR]; - amrex::Real lval = quadraticspline_1d(shapetypex, rx); + amrex::Real lval = quadraticspline_1d(periodic[0] ? 3 : shapetypex, rx); vals[0] += lval * nodaldata(i + l, 0, 0, comp0); vals[1] += lval * nodaldata(i + l, 0, 0, comp1); } @@ -918,7 +920,7 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], shapetypey = 3; amrex::Real ry = (xp[YDIR] - (plo[YDIR] + (j + m) * dx[YDIR])) / dx[YDIR]; - amrex::Real mval = quadraticspline_1d(shapetypey, ry); + amrex::Real mval = quadraticspline_1d(periodic[1] ? 3 : shapetypey, ry); for (int l = lmin; l < lmax; l++) { int shapetypex; @@ -932,7 +934,8 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], shapetypex = 3; amrex::Real rx = (xp[XDIR] - (plo[XDIR] + (i + l) * dx[XDIR])) / dx[XDIR]; - amrex::Real lval = quadraticspline_1d(shapetypex, rx); + amrex::Real lval = + quadraticspline_1d(periodic[0] ? 3 : shapetypex, rx); amrex::Real weight = lval * mval; vals[0] += weight * nodaldata(i + l, j + m, 0, comp0); vals[1] += weight * nodaldata(i + l, j + m, 0, comp1); @@ -962,7 +965,7 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], shapetypez = 3; amrex::Real rz = (xp[ZDIR] - (plo[ZDIR] + (k + n) * dx[ZDIR])) / dx[ZDIR]; - amrex::Real nval = quadraticspline_1d(shapetypez, rz); + amrex::Real nval = quadraticspline_1d(periodic[2] ? 3 : shapetypez, rz); for (int m = mmin; m < mmax; m++) { int shapetypey; @@ -976,7 +979,8 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], shapetypey = 3; amrex::Real ry = (xp[YDIR] - (plo[YDIR] + (j + m) * dx[YDIR])) / dx[YDIR]; - amrex::Real mval = quadraticspline_1d(shapetypey, ry); + amrex::Real mval = + quadraticspline_1d(periodic[1] ? 3 : shapetypey, ry); for (int l = lmin; l < lmax; l++) { int shapetypex; @@ -990,7 +994,8 @@ quadratic_interp_two(amrex::Real xp[AMREX_SPACEDIM], shapetypex = 3; amrex::Real rx = (xp[XDIR] - (plo[XDIR] + (i + l) * dx[XDIR])) / dx[XDIR]; - amrex::Real lval = quadraticspline_1d(shapetypex, rx); + amrex::Real lval = + quadraticspline_1d(periodic[0] ? 3 : shapetypex, rx); amrex::Real weight = lval * mval * nval; vals[0] += weight * nodaldata(i + l, j + m, k + n, comp0); vals[1] += weight * nodaldata(i + l, j + m, k + n, comp1); @@ -1029,6 +1034,7 @@ cubic_interp(amrex::Real xp[AMREX_SPACEDIM], const amrex::GpuArray dx, amrex::Array4 nodaldata, int comp, + const amrex::GpuArray &periodic, const int *lo, const int *hi) { @@ -1058,7 +1064,7 @@ cubic_interp(amrex::Real xp[AMREX_SPACEDIM], amrex::Real cell_center = plo[0] + (i + l) * dx[0]; amrex::Real rx = (xp[0] - cell_center) * inv_dx; - amrex::Real lval = cubicspline_1d(shapetypex, rx); + amrex::Real lval = cubicspline_1d(periodic[0] ? 3 : shapetypex, rx); value += lval * nodaldata(i + l, 0, 0, comp); } @@ -1107,8 +1113,8 @@ cubic_interp(amrex::Real xp[AMREX_SPACEDIM], amrex::Real rx = (xp[0] - (cell_centerx)) * inv_dx; amrex::Real ry = (xp[1] - (cell_centery)) * inv_dy; - amrex::Real lval = cubicspline_1d(shapetypex, rx); - amrex::Real mval = cubicspline_1d(shapetypey, ry); + amrex::Real lval = cubicspline_1d(periodic[0] ? 3 : shapetypex, rx); + amrex::Real mval = cubicspline_1d(periodic[1] ? 3 : shapetypey, ry); value += lval * mval * nodaldata(i + l, j + m, 0, comp); } @@ -1178,9 +1184,12 @@ cubic_interp(amrex::Real xp[AMREX_SPACEDIM], amrex::Real ry = (xp[1] - (cell_centery)) * inv_dy; amrex::Real rz = (xp[2] - (cell_centerz)) * inv_dz; - amrex::Real lval = cubicspline_1d(shapetypex, rx); - amrex::Real mval = cubicspline_1d(shapetypey, ry); - amrex::Real nval = cubicspline_1d(shapetypez, rz); + amrex::Real lval = + cubicspline_1d(periodic[0] ? 3 : shapetypex, rx); + amrex::Real mval = + cubicspline_1d(periodic[1] ? 3 : shapetypey, ry); + amrex::Real nval = + cubicspline_1d(periodic[2] ? 3 : shapetypez, rz); value += lval * mval * nval * nodaldata(i + l, j + m, k + n, comp); @@ -1223,6 +1232,7 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], amrex::Array4 nodaldata, int comp0, int comp1, + const amrex::GpuArray &periodic, const int *lo, const int *hi, amrex::Real vals[2]) @@ -1247,7 +1257,7 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], else shapetypex = 3; amrex::Real rx = (xp[0] - (plo[0] + (i + l) * dx[0])) * inv_dx; - amrex::Real lval = cubicspline_1d(shapetypex, rx); + amrex::Real lval = cubicspline_1d(periodic[0] ? 3 : shapetypex, rx); vals[0] += lval * nodaldata(i + l, 0, 0, comp0); vals[1] += lval * nodaldata(i + l, 0, 0, comp1); } @@ -1273,7 +1283,7 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], else shapetypey = 3; amrex::Real ry = (xp[1] - (plo[1] + (j + m) * dx[1])) * inv_dy; - amrex::Real mval = cubicspline_1d(shapetypey, ry); + amrex::Real mval = cubicspline_1d(periodic[1] ? 3 : shapetypey, ry); for (int l = lmin; l < lmax; ++l) { int shapetypex; @@ -1286,7 +1296,7 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], else shapetypex = 3; amrex::Real rx = (xp[0] - (plo[0] + (i + l) * dx[0])) * inv_dx; - amrex::Real lval = cubicspline_1d(shapetypex, rx); + amrex::Real lval = cubicspline_1d(periodic[0] ? 3 : shapetypex, rx); amrex::Real weight = lval * mval; vals[0] += weight * nodaldata(i + l, j + m, 0, comp0); vals[1] += weight * nodaldata(i + l, j + m, 0, comp1); @@ -1318,7 +1328,7 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], else shapetypez = 3; amrex::Real rz = (xp[2] - (plo[2] + (k + n) * dx[2])) * inv_dz; - amrex::Real nval = cubicspline_1d(shapetypez, rz); + amrex::Real nval = cubicspline_1d(periodic[2] ? 3 : shapetypez, rz); for (int m = mmin; m < mmax; ++m) { int shapetypey; @@ -1331,7 +1341,7 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], else shapetypey = 3; amrex::Real ry = (xp[1] - (plo[1] + (j + m) * dx[1])) * inv_dy; - amrex::Real mval = cubicspline_1d(shapetypey, ry); + amrex::Real mval = cubicspline_1d(periodic[1] ? 3 : shapetypey, ry); for (int l = lmin; l < lmax; ++l) { int shapetypex; @@ -1344,7 +1354,8 @@ cubic_interp_two(amrex::Real xp[AMREX_SPACEDIM], else shapetypex = 3; amrex::Real rx = (xp[0] - (plo[0] + (i + l) * dx[0])) * inv_dx; - amrex::Real lval = cubicspline_1d(shapetypex, rx); + amrex::Real lval = + cubicspline_1d(periodic[0] ? 3 : shapetypex, rx); amrex::Real weight = lval * mval * nval; vals[0] += weight * nodaldata(i + l, j + m, k + n, comp0); vals[1] += weight * nodaldata(i + l, j + m, k + n, comp1); diff --git a/Source/mpm_particle_grid_ops.cpp b/Source/mpm_particle_grid_ops.cpp index f931632..48bf20d 100644 --- a/Source/mpm_particle_grid_ops.cpp +++ b/Source/mpm_particle_grid_ops.cpp @@ -1012,17 +1012,19 @@ void MPMParticleContainer::interpolate_from_grid( } else if (order_scheme_directional[dim] == 2) { - quadratic_interp_two( - xp, iv, min_index, max_index, plo, dx, - nodal_data_arr, VELX_INDEX + dim, - DELTA_VELX_INDEX + dim, lo, hi, interp_vals); + quadratic_interp_two(xp, iv, min_index, max_index, + plo, dx, nodal_data_arr, + VELX_INDEX + dim, + DELTA_VELX_INDEX + dim, + periodic, lo, hi, interp_vals); } else if (order_scheme_directional[dim] == 3) { - cubic_interp_two( - xp, iv, min_index, max_index, plo, dx, - nodal_data_arr, VELX_INDEX + dim, - DELTA_VELX_INDEX + dim, lo, hi, interp_vals); + cubic_interp_two(xp, iv, min_index, max_index, plo, + dx, nodal_data_arr, + VELX_INDEX + dim, + DELTA_VELX_INDEX + dim, periodic, + lo, hi, interp_vals); } p.rdata(realData::xvel_prime + dim) = interp_vals[0]; @@ -1209,13 +1211,15 @@ void MPMParticleContainer::interpolate_from_grid_temperature( { p.rdata(realData::temperature) += quadratic_interp( xp, iv, min_index, max_index, plo, dx, - nodal_data_arr, DELTA_TEMPERATURE, lo, hi); + nodal_data_arr, DELTA_TEMPERATURE, periodic, lo, + hi); } else if (order_scheme_directional[0] == 3) { - p.rdata(realData::temperature) += cubic_interp( - xp, iv, min_index, max_index, plo, dx, - nodal_data_arr, DELTA_TEMPERATURE, lo, hi); + p.rdata(realData::temperature) += + cubic_interp(xp, iv, min_index, max_index, plo, dx, + nodal_data_arr, DELTA_TEMPERATURE, + periodic, lo, hi); } } From 28b3b35337972591b8d561c1e12af6613fbe73f0 Mon Sep 17 00:00:00 2001 From: SreejithNREL Date: Thu, 13 Aug 2026 13:18:07 -0600 Subject: [PATCH 06/11] bug fix for heatflux+convective temperature test cases --- .gitignore | 2 ++ Developer_Tools/Run_All_Tests.py | 22 +++++++++---------- Source/mpm_init.cpp | 1 - Source/nodal_data_ops.cpp | 3 +-- Source/utilities.cpp | 5 ++--- .../PreProcess/config.json | 4 ++-- .../2D_Heat_Conduction/PreProcess/config.json | 2 +- Tests/Dam_Break/PreProcess/config.json | 4 ++-- .../PreProcess/config.json | 1 + 9 files changed, 22 insertions(+), 22 deletions(-) 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 df6def7..6c8b776 100644 --- a/Developer_Tools/Run_All_Tests.py +++ b/Developer_Tools/Run_All_Tests.py @@ -2518,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], @@ -2546,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], @@ -2571,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], @@ -2620,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], @@ -2670,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], @@ -2843,25 +2843,25 @@ def _run_parameter_sweeps(): if test_name == "1D_Axial_Bar_Vibration": print('Nothing to do') - #Run_ParameterSweep_1D_Axial_Bar_Vibration(cfg) + 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) elif test_name == "Dam_Break": print('Nothing to do') - #Run_ParameterSweep_Dambreak(cfg) + Run_ParameterSweep_Dambreak(cfg) elif test_name == "Elastic_disk_collision": print('Nothing to do') Run_ParameterSweep_EDC(cfg) diff --git a/Source/mpm_init.cpp b/Source/mpm_init.cpp index db763f4..f9d4fed 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 diff --git a/Source/nodal_data_ops.cpp b/Source/nodal_data_ops.cpp index e0f6f97..b6b7416 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) { diff --git a/Source/utilities.cpp b/Source/utilities.cpp index 0026f98..6007763 100644 --- a/Source/utilities.cpp +++ b/Source/utilities.cpp @@ -247,8 +247,8 @@ void Apply_Nodal_BCs_Temperature(amrex::Geometry &geom, 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; + bclo_dirichlet[d] = specs.bclo_temp[d]; + bchi_dirichlet[d] = specs.bchi_temp[d]; } nodal_bcs_temperature( geom, nodaldata, bclo_dirichlet.data(), bchi_dirichlet.data(), @@ -520,7 +520,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_Heat_Conduction_HeatFlux/PreProcess/config.json b/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json index 090696e..bc84798 100644 --- a/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction_HeatFlux/PreProcess/config.json @@ -40,13 +40,13 @@ "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/config.json b/Tests/2D_Heat_Conduction/PreProcess/config.json index 775a880..0358d6e 100644 --- a/Tests/2D_Heat_Conduction/PreProcess/config.json +++ b/Tests/2D_Heat_Conduction/PreProcess/config.json @@ -40,7 +40,7 @@ "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": "hdf5", "materialpoint_filename": "mpm_particles.h5", diff --git a/Tests/Dam_Break/PreProcess/config.json b/Tests/Dam_Break/PreProcess/config.json index a50a190..5af81f8 100644 --- a/Tests/Dam_Break/PreProcess/config.json +++ b/Tests/Dam_Break/PreProcess/config.json @@ -46,7 +46,7 @@ "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", @@ -89,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/Elastic_disk_collision/PreProcess/config.json b/Tests/Elastic_disk_collision/PreProcess/config.json index 521196e..b04558e 100644 --- a/Tests/Elastic_disk_collision/PreProcess/config.json +++ b/Tests/Elastic_disk_collision/PreProcess/config.json @@ -94,6 +94,7 @@ 0.0, 0.0 ], + "density": 1.0, "boundary_conditions": { "xlo": { "mom": "periodic" From 8d3726a40b4ad463366c88bc28f9de21605a25e7 Mon Sep 17 00:00:00 2001 From: SreejithNREL Date: Thu, 13 Aug 2026 13:32:54 -0600 Subject: [PATCH 07/11] cleaned --- Source/mpm_specs.H | 12 +++++++++--- Source/nodal_data_ops.cpp | 2 +- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/Source/mpm_specs.H b/Source/mpm_specs.H index 2d13fae..1f7c7b9 100644 --- a/Source/mpm_specs.H +++ b/Source/mpm_specs.H @@ -517,13 +517,19 @@ class MPMspecs if (pp.query(face_lo_temp_keys[d].c_str(), bc_temp_str)) { bclo_temp[d] = bc_temp_str_to_int(bc_temp_str); - amrex::Print()<<"\n bc low string = "< Date: Thu, 13 Aug 2026 13:47:05 -0600 Subject: [PATCH 08/11] a few more bug fixes --- Source/mpm_specs.H | 9 --------- Source/nodal_data_ops.cpp | 6 ++++-- 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/Source/mpm_specs.H b/Source/mpm_specs.H index 1f7c7b9..8a44a46 100644 --- a/Source/mpm_specs.H +++ b/Source/mpm_specs.H @@ -517,19 +517,10 @@ class MPMspecs if (pp.query(face_lo_temp_keys[d].c_str(), bc_temp_str)) { bclo_temp[d] = bc_temp_str_to_int(bc_temp_str); - amrex::Print() - << "\n bc low string = " << face_lo_temp_keys[d].c_str() - << " " << bc_temp_str << " " << bclo_temp[d] << " " << d - << "\n"; } if (pp.query(face_hi_temp_keys[d].c_str(), bc_temp_str)) { bchi_temp[d] = bc_temp_str_to_int(bc_temp_str); - amrex::Print() - << "\n bc hi string = " << face_hi_temp_keys[d].c_str() - << " " << bc_temp_str << " " << bchi_temp[d] << " " - << BC_TEMP_CONVECTION << " " << d << "\n"; - // amrex::Abort("\n Force exiting"); } ParmParse pp_lo_temp(face_lo_temp_ns[d]); diff --git a/Source/nodal_data_ops.cpp b/Source/nodal_data_ops.cpp index eb45e6e..b405ecf 100644 --- a/Source/nodal_data_ops.cpp +++ b/Source/nodal_data_ops.cpp @@ -895,7 +895,8 @@ void nodal_bcs_temperature(const amrex::Geometry geom, amrex::Real hc = is_lo ? h_lo_g[d] : h_hi_g[d]; 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 = h_conv_v * dx[dom_dir] / k_node; arr(nodeid, TEMPERATURE) = (arr(nb, TEMPERATURE) + Bi * Tinf) / (1.0 + Bi); bc_applied = true; @@ -1430,7 +1431,8 @@ 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 k_node = (m > shunya) ? mk / m : eka; + amrex::Real Bi = h_conv_v * dx[dom_dir] / k_node; arr(nodeid, TEMPERATURE) = (arr(nb, TEMPERATURE) + Bi * Tinf) / (eka + Bi); } From f9a309c413bb1d8d23498ff3c008a3652e47ebce Mon Sep 17 00:00:00 2001 From: Sreejith Date: Thu, 13 Aug 2026 14:25:55 -0600 Subject: [PATCH 09/11] a few few more bug fixes --- Developer_Tools/Run_All_Tests.py | 8 +++--- Source/nodal_data_ops.cpp | 8 ++++-- Source/utilities.cpp | 47 +++++++------------------------- 3 files changed, 20 insertions(+), 43 deletions(-) diff --git a/Developer_Tools/Run_All_Tests.py b/Developer_Tools/Run_All_Tests.py index 6c8b776..bbdd3ef 100644 --- a/Developer_Tools/Run_All_Tests.py +++ b/Developer_Tools/Run_All_Tests.py @@ -2843,7 +2843,7 @@ def _run_parameter_sweeps(): if test_name == "1D_Axial_Bar_Vibration": print('Nothing to do') - Run_ParameterSweep_1D_Axial_Bar_Vibration(cfg) + #Run_ParameterSweep_1D_Axial_Bar_Vibration(cfg) elif test_name == "1D_Heat_Conduction": print('Nothing to do') Run_ParameterSweep_1D_HeatConduction(cfg) @@ -2858,13 +2858,13 @@ def _run_parameter_sweeps(): Run_ParameterSweep_2D_HeatConduction(cfg) elif test_name == "2D_Heat_Conduction_Cylinder_Dirichlet": print('Nothing to do') - #Run_ParameterSweep_2D_HeatConduction_Cylinder_Dirichlet(cfg) + Run_ParameterSweep_2D_HeatConduction_Cylinder_Dirichlet(cfg) elif test_name == "Dam_Break": print('Nothing to do') - Run_ParameterSweep_Dambreak(cfg) + #Run_ParameterSweep_Dambreak(cfg) elif test_name == "Elastic_disk_collision": print('Nothing to do') - Run_ParameterSweep_EDC(cfg) + #Run_ParameterSweep_EDC(cfg) # Save results _sweep_results_path = os.path.join(ROOT, "sweep_results.json") diff --git a/Source/nodal_data_ops.cpp b/Source/nodal_data_ops.cpp index b405ecf..a0765a8 100644 --- a/Source/nodal_data_ops.cpp +++ b/Source/nodal_data_ops.cpp @@ -893,10 +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 k_node = (m > shunya) ? mk / m : eka; - amrex::Real Bi = h_conv_v * dx[dom_dir] / k_node; + amrex::Real Bi = hc * dx_g[d] / k_node; arr(nodeid, TEMPERATURE) = (arr(nb, TEMPERATURE) + Bi * Tinf) / (1.0 + Bi); bc_applied = true; @@ -1431,8 +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 mk = arr(nodeid, MASS_CONDUCTIVITY); + amrex::Real m = arr(nodeid, MASS_INDEX); amrex::Real k_node = (m > shunya) ? mk / m : eka; - amrex::Real Bi = h_conv_v * dx[dom_dir] / k_node; + 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 6007763..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]; - bchi_dirichlet[d] = specs.bchi_temp[d]; - } - 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 From da1bcd6312d0c93e6e57009f63f8ab38b1a20260 Mon Sep 17 00:00:00 2001 From: Sreejith Date: Thu, 13 Aug 2026 17:07:37 -0600 Subject: [PATCH 10/11] bug fixes from copilot review --- Source/mpm_init.cpp | 13 +++++++++++++ Source/mpm_particle_container.cpp | 7 ++++++- Source/mpm_particle_timestep.cpp | 11 ++++++++--- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- .../PreProcess/Generate_MPs_Inputfile_Generic.py | 8 +++++++- 12 files changed, 90 insertions(+), 13 deletions(-) diff --git a/Source/mpm_init.cpp b/Source/mpm_init.cpp index f9d4fed..b87be18 100644 --- a/Source/mpm_init.cpp +++ b/Source/mpm_init.cpp @@ -812,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]; diff --git a/Source/mpm_particle_container.cpp b/Source/mpm_particle_container.cpp index 99082fb..9b9ecee 100644 --- a/Source/mpm_particle_container.cpp +++ b/Source/mpm_particle_container.cpp @@ -202,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 @@ -255,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 b36b192..c95c72d 100644 --- a/Source/mpm_particle_timestep.cpp +++ b/Source/mpm_particle_timestep.cpp @@ -55,13 +55,13 @@ amrex::Real MPMParticleContainer::Calculate_time_step(MPMspecs &specs) if (p.idata(intData::phase) == 0) { amrex::Real Cs = 0.0; - if (p.idata(intData::constitutive_model) == 1 or - p.idata(intData::constitutive_model) == 2) + if (p.idata(intData::constitutive_model) == 1) { 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) * @@ -73,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/Tests/1D_Axial_Bar_Vibration/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/1D_Axial_Bar_Vibration/PreProcess/Generate_MPs_Inputfile_Generic.py index e480a38..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 @@ -260,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"} 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 e480a38..e75f523 100644 --- a/Tests/1D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/1D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -260,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"} 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 e480a38..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 @@ -260,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"} 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 e480a38..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 @@ -260,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"} 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 e480a38..e75f523 100644 --- a/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/2D_Heat_Conduction/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -260,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"} 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 e480a38..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 @@ -260,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"} diff --git a/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py b/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py index e480a38..5341d92 100644 --- a/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/Dam_Break/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -260,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"} 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 e480a38..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 @@ -260,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"} 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 e480a38..575c69d 100644 --- a/Tests/Elastic_disk_collision/PreProcess/Generate_MPs_Inputfile_Generic.py +++ b/Tests/Elastic_disk_collision/PreProcess/Generate_MPs_Inputfile_Generic.py @@ -260,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"} From c12be0ab51b8dfff5f3b70cafde7dba703bfaabf Mon Sep 17 00:00:00 2001 From: Sreejith Date: Thu, 13 Aug 2026 17:23:14 -0600 Subject: [PATCH 11/11] changed back to ascii for test cases in ci --- Tests/1D_Axial_Bar_Vibration/PreProcess/config.json | 6 +++--- Tests/1D_Heat_Conduction/PreProcess/config.json | 6 +++--- .../1D_Heat_Conduction_Convective/PreProcess/config.json | 8 ++++---- Tests/2D_Heat_Conduction/PreProcess/config.json | 6 +++--- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json b/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json index e2261f8..919137c 100644 --- a/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json +++ b/Tests/1D_Axial_Bar_Vibration/PreProcess/config.json @@ -45,9 +45,9 @@ "stress_update_scheme": "MUSL", "output_tag": "1D_Axial_Bar_Vibration_dim1_npcx2_ord1_flip1.0_susMUSL_CFL0.1_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPFalse_BuildSystemcmake_44edc4", "input_filename": "Inputs_1DAxialBarVibration.inp", - "build_with_hdf": true, - "output_format": "hdf5", - "materialpoint_filename": "mpm_particles.h5", + "build_with_hdf": false, + "output_format": "ascii", + "materialpoint_filename": "mpm_particles.dat", "plot_to_check": true, "build_system": "cmake", "use_mpi": true, diff --git a/Tests/1D_Heat_Conduction/PreProcess/config.json b/Tests/1D_Heat_Conduction/PreProcess/config.json index cde3982..47a1149 100644 --- a/Tests/1D_Heat_Conduction/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction/PreProcess/config.json @@ -42,10 +42,10 @@ "stress_update_scheme": "MUSL", "output_tag": "1D_Heat_Conduction_npcx2_ord1_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemcmake_868edf", "input_filename": "Inputs_1DHeatConduction.inp", - "output_format": "hdf5", - "materialpoint_filename": "mpm_particles.h5", + "output_format": "ascii", + "materialpoint_filename": "mpm_particles.dat", "plot_to_check": true, - "build_with_hdf": true, + "build_with_hdf": false, "build_system": "cmake", "use_mpi": true, "use_cuda": false, diff --git a/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json b/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json index 7fdca15..c8064a3 100644 --- a/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json +++ b/Tests/1D_Heat_Conduction_Convective/PreProcess/config.json @@ -9,8 +9,8 @@ "ny": 4 }, "ppc": [ - 1, - 1 + 2, + 2 ], "bodies": [ { @@ -38,7 +38,7 @@ "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_ord2_susMUSL_USEHDFTrue_OFORMhdf5_MPI=True_CUDAFalse_HIPFalse_OMPFalse_SYCLFalse_EBFalse_TEMPTrue_BuildSystemgnumake_9af11f", "input_filename": "Inputs_1DHeatConduction_Convective.inp", @@ -105,4 +105,4 @@ "do_calculate_minmaxpos": 0, "write_diag_output_time": 0.01 } -} +} \ No newline at end of file diff --git a/Tests/2D_Heat_Conduction/PreProcess/config.json b/Tests/2D_Heat_Conduction/PreProcess/config.json index 0358d6e..3363f1d 100644 --- a/Tests/2D_Heat_Conduction/PreProcess/config.json +++ b/Tests/2D_Heat_Conduction/PreProcess/config.json @@ -42,10 +42,10 @@ "stress_update_scheme": "MUSL", "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": "hdf5", - "materialpoint_filename": "mpm_particles.h5", + "output_format": "ascii", + "materialpoint_filename": "mpm_particles.dat", "plot_to_check": false, - "build_with_hdf": true, + "build_with_hdf": false, "build_system": "cmake", "use_mpi": true, "use_cuda": false,