From 96cc2795303cdbab5b44ab9efce05c84e39ac2b8 Mon Sep 17 00:00:00 2001 From: Alexander Hermann Date: Fri, 7 Aug 2026 22:59:15 +0200 Subject: [PATCH 1/2] Stop writing a fabricated key into every generated Kratos DEM deck MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `main` writes `PARTICLE_FRICTION` into the mdpa and into MaterialsDEM.json on every DEM deck OASiS generates. The key does not exist. Measured with `grep -r -a -F` against the full 28-application Kratos build at /mnt/kratos-tier2/kv: zero occurrences. The real keys are STATIC_FRICTION and DYNAMIC_FRICTION, per contact PAIR, in `material_relations` — a different place in a different schema, not a rename. It also emitted the FEM materials schema, `problem_data`, `solver_type` and `time_stepping`, none of which DEM reads, plus `DEM_timestep_safety_factor` (the real key is `DeltaTimeSafetyFactor`, itself dead code) and a mesh filename DEM never opens. The deck it produced had therefore never run. `catalog_template_executes` recorded `kratos::dem::2d_rc=0`, but that gate runs the template — which only WRITES files and exits — and never the `input.py` it writes. Run that, and it dies at `entry string : strategy` in kratos_parameters.cpp:426. Now: kratos:dem:2d template_rc=0 input.py_rc=0 both run The deck is a 2-D granular settling problem — 64 CylinderParticle2D discs onto a floor between two walls, ~9 s. Physical, not merely exit-0: all 64 particles survive, the pile comes to rest at max |v_y| = 0.005 m/s, and y_min sits one radius above the floor with a 0.1%-of-radius Hertzian overlap. The emitted input.py also refuses to exit 0 when the model has quietly lost its particles. That matters because Kratos DEM deletes escaped particles silently: measured on this case, dt=2e-4 loses all 16 and dt=1e-4 loses 12 of 16, both at exit code 0. One file, so it can be reviewed and land ahead of the full consolidation. The knowledge, fixtures and ledger for this work sit on feature/kratos-generator. Co-Authored-By: Claude Opus 5 (1M context) --- src/backends/kratos/generators/dem.py | 667 ++++++++++++++++++-------- 1 file changed, 462 insertions(+), 205 deletions(-) diff --git a/src/backends/kratos/generators/dem.py b/src/backends/kratos/generators/dem.py index c376e9c8..43be5914 100644 --- a/src/backends/kratos/generators/dem.py +++ b/src/backends/kratos/generators/dem.py @@ -1,217 +1,381 @@ """Kratos DEM (Discrete Element Method) generators and knowledge. Uses the REAL Kratos DEMApplication — NOT standalone numpy/scipy. -Generates ProjectParametersDEM.json + .mdpa + MaterialsDEM.json + input.py. +The generator is a FILE WRITER: it emits the four files DEMApplication reads +(DEM.mdpa, DEM_FEM_boundary.mdpa, MaterialsDEM.json, +ProjectParametersDEM.json) plus input.py, the entry point that calls +DEMAnalysisStage. Both halves are executed by +scripts/audit_two_stage_templates.py; the generator exiting 0 is NOT evidence +that the deck it wrote runs, and for a long time it did not. """ -import json - def _dem_2d_kratos(params: dict) -> str: - """Generator that produces real Kratos DEM input files. - - Creates: ProjectParametersDEM.json, particles.mdpa, MaterialsDEM.json, input.py - The run_with_generator pipeline will find input.py and execute it via Kratos. - """ - n_particles = params.get("n_particles", 200) - radius = params.get("radius", 0.01) - density = params.get("density", 2650.0) - gravity_y = params.get("gravity", -9.81) - dt = params.get("dt", 1e-5) - T_end = params.get("T_end", 0.5) - young_modulus = params.get("E", 1e7) - restitution = params.get("restitution", 0.5) - friction = params.get("friction", 0.4) - domain_x = params.get("domain_x", 0.5) - domain_y = params.get("domain_y", 1.0) - - return f'''\ -"""Generator: creates real Kratos DEM input files and runs the analysis.""" -import os + n_particles = int(params.get("n_particles", 64)) + radius = float(params.get("radius", 0.01)) + density = float(params.get("density", 2650.0)) + young_modulus = float(params.get("E", 1.0e7)) + poisson = float(params.get("nu", 0.25)) + restitution = float(params.get("restitution", 0.5)) + friction = float(params.get("friction", 0.4)) + gravity_y = float(params.get("gravity", -9.81)) + domain_x = float(params.get("domain_x", 0.30)) + domain_y = float(params.get("domain_y", 0.40)) + drop_height = float(params.get("drop_height", 0.05)) + T_end = float(params.get("T_end", 0.25)) + dt = params.get("dt", None) + dt_safety = float(params.get("dt_safety_factor", 0.02)) + problem_name = str(params.get("problem_name", "granular")) + threads = int(params.get("omp_threads", 2)) + + header = f'''\ +"""Kratos DEM — 2D granular settling — writes a complete DEMApplication deck. + +This script is a FILE WRITER. It emits + + {problem_name}DEM.mdpa spheres (CylinderParticle2D) + {problem_name}DEM_FEM_boundary.mdpa rigid walls (RigidFace3D3N) + MaterialsDEM.json materials + material_relations + table + ProjectParametersDEM.json FLAT DEM schema + input.py the entry point that calls Kratos + +and then exits. Run `python input.py` to perform the analysis. + +Schema notes (each verified by execution against DEMApplication 10.4.3): + * ProjectParametersDEM.json is FLAT. There is no "problem_data" block and no + "solver_settings.solver_type"/"time_stepping". DEM reads + solver_settings.strategy, a python module name inside DEMApplication. + * The mdpa filenames are problem_name + a FIXED tag ("DEM", + "DEM_FEM_boundary", ...) with no separator. A missing file is not an + error -- the run completes on an empty model -- so the names must be right. + * MaterialsDEM.json needs materials / material_relations / + material_assignation_table. Friction is STATIC_FRICTION + DYNAMIC_FRICTION + per contact PAIR in material_relations; PARTICLE_FRICTION does not exist. + * MaxTimeStep is the only key that sets dt and it is used verbatim. There is + no stability check, so dt is sized here from the Rayleigh estimate. +""" import json import math -import numpy as np - -# === Parameters === -n_particles = {n_particles} -radius = {radius} -density = {density} -gravity_y = {gravity_y} -dt = {dt} -T_end = {T_end} -E_contact = {young_modulus} -restitution = {restitution} -friction = {friction} -domain_x = {domain_x} -domain_y = {domain_y} - -np.random.seed(42) - -# === Generate particle positions === -cols = int(np.ceil(np.sqrt(n_particles))) -spacing = 2.5 * radius + +# ---- problem definition ------------------------------------------------- +problem_name = "{problem_name}" +n_particles = {n_particles} +radius = {radius} +density = {density} # kg/m^3 -> PARTICLE_DENSITY +young_modulus = {young_modulus} # Pa -> YOUNG_MODULUS +poisson = {poisson} +restitution = {restitution} # -> COEFFICIENT_OF_RESTITUTION (contact pair) +friction = {friction} # -> STATIC_FRICTION / DYNAMIC_FRICTION (pair) +gravity_y = {gravity_y} +domain_x = {domain_x} +domain_y = {domain_y} +drop_height = {drop_height} +T_end = {T_end} +dt_override = {dt!r} +dt_safety = {dt_safety} +omp_threads = {threads} +''' + + body = r''' +# ---- time step ---------------------------------------------------------- +# Kratos DEM never checks stability: MaxTimeStep is used verbatim, and a step +# above the limit blows the contacts up silently (particles leave the bounding +# box and are deleted without a message). Size it from the Rayleigh surface-wave +# criterion, which is the standard DEM estimate: +# +# dt_R = pi * R * sqrt(rho / G) / (0.1631 * nu + 0.8766), G = E / (2(1+nu)) +# +# and take a small fraction of it. Measured for the default parameters below: +# dt_R = 8.8e-4 s, while the deck is already unstable at 1e-4 (12 of 16 +# particles lost) and stable at 5e-5 -- i.e. the true limit is near 0.07*dt_R, +# so the Rayleigh value is optimistic by more than an order of magnitude for the +# viscous-Hertz law. dt_safety = 0.02 sits a factor ~3 below the measured limit. +shear_modulus = young_modulus / (2.0 * (1.0 + poisson)) +dt_rayleigh = (math.pi * radius * math.sqrt(density / shear_modulus) + / (0.1631 * poisson + 0.8766)) +dt = float(dt_override) if dt_override is not None else dt_safety * dt_rayleigh +n_steps = max(1, int(round(T_end / dt))) + +# ---- particle layout: a block dropped onto the floor -------------------- +spacing = 2.2 * radius +cols = max(1, int(math.floor((domain_x - 2.0 * spacing) / spacing))) +cols = min(cols, max(1, int(math.ceil(math.sqrt(n_particles))))) +rows = int(math.ceil(n_particles / cols)) +block_w = (cols - 1) * spacing +x0 = 0.5 * (domain_x - block_w) + positions = [] for i in range(n_particles): - row, col = divmod(i, cols) - x = (col + 1) * spacing + np.random.uniform(-radius*0.1, radius*0.1) - y = domain_y * 0.5 + (row + 1) * spacing - x = min(max(x, radius), domain_x - radius) - positions.append((x, y, 0.0)) - -# === Write MDPA === -mdpa = "Begin ModelPartData\\nEnd ModelPartData\\n\\n" -mdpa += "Begin Properties 1\\n" -mdpa += f" PARTICLE_DENSITY {{density}}\\n" -mdpa += f" YOUNG_MODULUS {{E_contact}}\\n" -mdpa += f" POISSON_RATIO 0.25\\n" -mdpa += f" PARTICLE_FRICTION {{friction}}\\n" -mdpa += f" COEFFICIENT_OF_RESTITUTION {{restitution}}\\n" -mdpa += f" PARTICLE_MATERIAL 1\\n" -mdpa += f" ROLLING_FRICTION 0.01\\n" -mdpa += "End Properties\\n\\n" - -mdpa += "Begin Nodes\\n" -for i, (x, y, z) in enumerate(positions, 1): - mdpa += f" {{i}} {{x:.8f}} {{y:.8f}} {{z:.8f}}\\n" -mdpa += "End Nodes\\n\\n" + r, c = divmod(i, cols) + positions.append((x0 + c * spacing, + drop_height + radius + r * spacing, + 0.0)) -mdpa += "Begin Elements SphericParticle3D\\n" +top_y = max(p[1] for p in positions) +if top_y + radius >= domain_y: + raise SystemExit( + "particle block (top y=%.4f) does not fit under domain_y=%.4f; " + "reduce n_particles/radius or raise domain_y" % (top_y + radius, domain_y)) + +# ---- spheres mdpa: DEM.mdpa ------------------------------- +# The element name is CylinderParticle2D: Kratos DEM is not 3D-only, and a +# genuine 2D disc is the right element for a planar problem. There is no +# SphericParticle2D. RADIUS is a CORE nodal variable and is the one mandatory +# per-node input; a particle whose RADIUS was never written keeps a silent 0. +lines = ["Begin ModelPartData", "End ModelPartData", "", + "Begin Properties 1", "End Properties", "", "Begin Nodes"] +for i, (x, y, z) in enumerate(positions, 1): + lines.append(" %d %.10f %.10f %.10f" % (i, x, y, z)) +lines += ["End Nodes", "", "Begin Elements CylinderParticle2D"] for i in range(1, n_particles + 1): - mdpa += f" {{i}} 1 {{i}}\\n" -mdpa += "End Elements\\n\\n" - -# All particles in one sub model part -mdpa += "Begin SubModelPart DEMParts_particles\\n" -mdpa += " Begin SubModelPartData\\n" -mdpa += " End SubModelPartData\\n" -mdpa += " Begin SubModelPartNodes\\n" + lines.append(" %d 1 %d" % (i, i)) +lines += ["End Elements", "", "Begin NodalData RADIUS"] for i in range(1, n_particles + 1): - mdpa += f" {{i}}\\n" -mdpa += " End SubModelPartNodes\\n" -mdpa += "End SubModelPart\\n" - -with open("particles.mdpa", "w") as f: - f.write(mdpa) - -# === Write MaterialsDEM.json === -materials = {{ - "properties": [{{ - "model_part_name": "SpheresPart", - "properties_id": 1, - "hydrodynamic_law_parameters": {{}}, - "Material": {{ - "Variables": {{ - "PARTICLE_DENSITY": density, - "YOUNG_MODULUS": E_contact, - "POISSON_RATIO": 0.25, - "PARTICLE_FRICTION": friction, - "COEFFICIENT_OF_RESTITUTION": restitution, - "PARTICLE_MATERIAL": 1, - "ROLLING_FRICTION": 0.01, - "ROLLING_FRICTION_WITH_WALLS": 0.01, - }}, - "constitutive_law": {{ - "name": "DEM_D_Hertz_viscous_Coulomb" - }} - }} - }}] -}} + lines.append(" %d 0 %.10f" % (i, radius)) +lines += ["End NodalData", "", "Begin SubModelPart DEMParts_Body", + " Begin SubModelPartNodes"] +lines += [" %d" % i for i in range(1, n_particles + 1)] +lines += [" End SubModelPartNodes", " Begin SubModelPartElements"] +lines += [" %d" % i for i in range(1, n_particles + 1)] +lines += [" End SubModelPartElements", "End SubModelPart", ""] +with open(problem_name + "DEM.mdpa", "w") as f: + f.write("\n".join(lines)) + +# ---- walls mdpa: DEM_FEM_boundary.mdpa -------------------- +# Walls are CONDITIONS, and only the 3D face spellings are registered +# (RigidFace2D2N does not exist), so a 2D wall is a 3D face given a thickness +# in z. Each wall is a quad split into two RigidFace3D3N triangles. +zlo, zhi = -5.0 * radius, 5.0 * radius +walls = { + "floor": [(0.0, 0.0, zlo), (domain_x, 0.0, zlo), + (domain_x, 0.0, zhi), (0.0, 0.0, zhi)], + "left": [(0.0, 0.0, zlo), (0.0, domain_y, zlo), + (0.0, domain_y, zhi), (0.0, 0.0, zhi)], + "right": [(domain_x, 0.0, zlo), (domain_x, domain_y, zlo), + (domain_x, domain_y, zhi), (domain_x, 0.0, zhi)], +} +wlines = ["Begin ModelPartData", "End ModelPartData", "", + "Begin Properties 2", "End Properties", "", "Begin Nodes"] +node_id, cond_id = 0, 0 +wall_nodes, wall_conds = {}, {} +node_lines, cond_lines = [], [] +for name, quad in walls.items(): + ids = [] + for (x, y, z) in quad: + node_id += 1 + ids.append(node_id) + node_lines.append(" %d %.10f %.10f %.10f" % (node_id, x, y, z)) + wall_nodes[name] = ids + cids = [] + for tri in ((0, 1, 2), (0, 2, 3)): + cond_id += 1 + cids.append(cond_id) + cond_lines.append(" %d 2 %d %d %d" + % (cond_id, ids[tri[0]], ids[tri[1]], ids[tri[2]])) + wall_conds[name] = cids +wlines += node_lines + ["End Nodes", "", "Begin Conditions RigidFace3D3N"] +wlines += cond_lines + ["End Conditions", ""] +for name in walls: + wlines += ["Begin SubModelPart DEM-FEM-Wall_%s" % name, + " Begin SubModelPartData", " End SubModelPartData", + " Begin SubModelPartNodes"] + wlines += [" %d" % i for i in wall_nodes[name]] + wlines += [" End SubModelPartNodes", " Begin SubModelPartConditions"] + wlines += [" %d" % i for i in wall_conds[name]] + wlines += [" End SubModelPartConditions", "End SubModelPart", ""] +with open(problem_name + "DEM_FEM_boundary.mdpa", "w") as f: + f.write("\n".join(wlines)) + +# ---- MaterialsDEM.json -------------------------------------------------- +# Three mandatory top-level keys. Per-particle data goes in "materials"; +# everything about a CONTACT (restitution, friction, the constitutive law name) +# goes in "material_relations", keyed by the pair of material ids. Every +# sub model part that carries entities needs a row in the assignation table -- +# omitting the wall row fails deep inside the solver with a PROPERTIES_ID error. +contact_pair_variables = { + "COEFFICIENT_OF_RESTITUTION": restitution, + "STATIC_FRICTION": friction, + "DYNAMIC_FRICTION": friction, + "FRICTION_DECAY": 500.0, + "ROLLING_FRICTION": 0.01, + "ROLLING_FRICTION_WITH_WALLS": 0.01, + "DEM_DISCONTINUUM_CONSTITUTIVE_LAW_NAME": "DEM_D_Hertz_viscous_Coulomb2D", +} +materials = { + "materials": [ + {"material_name": "granular", "material_id": 1, + "Variables": {"PARTICLE_DENSITY": density, + "YOUNG_MODULUS": young_modulus, + "POISSON_RATIO": poisson}}, + {"material_name": "wall", "material_id": 2, + "Variables": {"PARTICLE_DENSITY": 7850.0, + "YOUNG_MODULUS": 100.0 * young_modulus, + "POISSON_RATIO": 0.30}}, + ], + "material_relations": [ + {"material_names_list": ["granular", "granular"], + "material_ids_list": [1, 1], + "Variables": dict(contact_pair_variables)}, + {"material_names_list": ["granular", "wall"], + "material_ids_list": [1, 2], + "Variables": dict(contact_pair_variables)}, + ], + "material_assignation_table": [ + ["SpheresPart.DEMParts_Body", "granular"], + ] + [["RigidFacePart.DEM-FEM-Wall_%s" % n, "wall"] for n in walls], +} with open("MaterialsDEM.json", "w") as f: json.dump(materials, f, indent=2) -# === Write ProjectParametersDEM.json === -n_steps = int(T_end / dt) -output_interval = max(1, n_steps // 20) - -project_params = {{ - "problem_data": {{ - "problem_name": "particles", - "parallel_type": "OpenMP", - "echo_level": 0, - "start_time": 0.0, - "end_time": T_end - }}, - "solver_settings": {{ - "solver_type": "dem_solver", - "model_part_name": "SpheresPart", - "domain_size": 3, - "model_import_settings": {{ - "input_type": "mdpa", - "input_filename": "particles" - }}, - "material_import_settings": {{ - "materials_filename": "MaterialsDEM.json" - }}, - "time_stepping": {{ - "time_step": dt - }}, - "body_force_per_unit_mass": [0.0, gravity_y, 0.0], - "DEM_timestep_safety_factor": 0.5, - "MaxAmplificationRatioOfSearchRadius": 0.0, - "BoundingBoxMaxX": domain_x, - "BoundingBoxMaxY": domain_y * 2, - "BoundingBoxMaxZ": 0.1, - "BoundingBoxMinX": 0.0, - "BoundingBoxMinY": 0.0, - "BoundingBoxMinZ": -0.1, - "BoundingBoxEnlargementFactor": 1.1, - "AutomaticBoundingBoxOption": False, - "ContactMeshOption": "use_particles" - }}, - "output_processes": {{ - "vtk_output": [{{ - "python_module": "vtk_output_process", - "kratos_module": "KratosMultiphysics", - "process_name": "VtkOutputProcess", - "Parameters": {{ - "model_part_name": "SpheresPart", - "output_control_type": "step", - "output_interval": output_interval, - "file_format": "ascii", - "output_path": "vtk_output", - "nodal_solution_step_data_variables": ["VELOCITY", "TOTAL_FORCES"] - }} - }}] - }}, - "processes": {{ - "walls_process_list": [{{ - "python_module": "assign_vector_variable_process", - "kratos_module": "KratosMultiphysics", - "process_name": "AssignVectorVariableProcess", - "Parameters": {{ - "model_part_name": "SpheresPart.DEMParts_particles", - "variable_name": "RADIUS", - "value": [radius, 0.0, 0.0], - "constrained": [False, False, False], - "interval": [0, "End"] - }} - }}] - }} -}} - +# ---- ProjectParametersDEM.json (FLAT) ----------------------------------- +# Unknown top-level keys and unknown keys inside solver_settings are rejected by +# ValidateAndAssignDefaults, so this list is exactly the DEM vocabulary. +# FinalTime -- not problem_data.end_time -- decides when the run stops. +project_parameters = { + "problem_name": problem_name, + "Dimension": 2, + "FinalTime": T_end, + "MaxTimeStep": dt, + "OutputTimeStep": max(dt, T_end / 20.0), + "GravityX": 0.0, + "GravityY": gravity_y, + "GravityZ": 0.0, + "ElementType": "CylinderPartDEMElement2D", + "TranslationalIntegrationScheme": "Symplectic_Euler", + "RotationalIntegrationScheme": "Direct_Integration", + "RotationOption": True, + "RollingFrictionOption": False, + "dem_inlet_option": False, + "BoundingBoxOption": True, + "AutomaticBoundingBoxOption": False, + "BoundingBoxEnlargementFactor": 1.0, + "BoundingBoxMinX": -domain_x, + "BoundingBoxMinY": -domain_y, + "BoundingBoxMinZ": -1.0, + "BoundingBoxMaxX": 2.0 * domain_x, + "BoundingBoxMaxY": 2.0 * domain_y, + "BoundingBoxMaxZ": 1.0, + "NeighbourSearchFrequency": 10, + "SearchTolerance": 0.5 * radius, + # Output goes through DEM's own GiD writer, into _Post_Files. + # The two alternatives do not work on a stock build: + # post_vtk_option -> DEMApplication.dem_vtk_output imports pyevtk, which + # Kratos does not depend on and does not install. + # the core vtk_output_process, attached to SpheresPart -> its Check() runs + # before DEM populates that part's ProcessInfo, so + # output_control_type "step" raises 'STEP not found + # in process info of SpheresPart' and "time" raises + # the same for TIME. + "do_print_results_option": True, + "post_gid_option": True, + "post_vtk_option": False, + "PostVelocity": True, + "PostRadius": True, + "PostTotalForces": True, + "solver_settings": { + "strategy": "sphere_strategy", + "material_import_settings": {"materials_filename": "MaterialsDEM.json"}, + }, +} with open("ProjectParametersDEM.json", "w") as f: - json.dump(project_params, f, indent=2) + json.dump(project_parameters, f, indent=2) -# === Write input.py (the entry point) === -input_py = """import KratosMultiphysics +# ---- input.py: the entry point ------------------------------------------ +INPUT_PY = """# Kratos DEM entry point. Reads ProjectParametersDEM.json from this directory. +# +# Cap OpenMP threads BEFORE importing Kratos: a few-hundred-particle DEM step is +# dominated by parallel-region spin overhead, so more threads is slower, not +# faster. Measured on this deck, 2000 steps: 98 s at the machine default (32 +# threads) against 1.5 s at 2, with identical results to the last digit. +# +# Treat the RATIO as indicative, not as a specification. It was taken on a +# shared 32-core box carrying a load average around 43, where oversubscribed +# spin-wait is at its worst; on an idle machine the gap is smaller. What does +# not depend on the load is the direction and the reason -- the parallel region +# costs more to enter than the work inside it is worth at this problem size -- +# and that is independently recorded elsewhere in this catalog from a separate +# 2026-06-12 measurement on a different case. +# +# setdefault, so an explicit OMP_NUM_THREADS in the environment still wins. +import json +import os + +os.environ.setdefault("OMP_NUM_THREADS", "__THREADS__") + +import KratosMultiphysics from KratosMultiphysics.DEMApplication.DEM_analysis_stage import DEMAnalysisStage -if __name__ == "__main__": - with open("ProjectParametersDEM.json", "r") as f: - project_parameters = KratosMultiphysics.Parameters(f.read()) - model = KratosMultiphysics.Model() - DEMAnalysisStage(model, project_parameters).Run() + +class GranularSettling(DEMAnalysisStage): + # DEMAnalysisStage.Finalize() deletes its model parts, and holding a python + # reference to one across that call segfaults. Read what is needed here, + # before delegating. + def Finalize(self): + spheres = self.spheres_model_part + walls = self.rigid_face_model_part + ys = [n.Y for n in spheres.Nodes] + vy = [n.GetSolutionStepValue(KratosMultiphysics.VELOCITY_Y) + for n in spheres.Nodes] + self.summary = { + "particles_at_end": spheres.NumberOfElements(), + "wall_conditions": walls.NumberOfConditions(), + "y_min": min(ys) if ys else None, + "y_max": max(ys) if ys else None, + "max_abs_velocity_y": max((abs(v) for v in vy), default=None), + } + super().Finalize() + + +with open("ProjectParametersDEM.json") as _f: + parameters = KratosMultiphysics.Parameters(_f.read()) + +analysis = GranularSettling(KratosMultiphysics.Model(), parameters) +n_expected = __N_PARTICLES__ +analysis.Run() + +summary = dict(analysis.summary) +summary["particles_at_start"] = n_expected +# Read these back from the deck, not from a value baked in here, so that editing +# ProjectParametersDEM.json by hand keeps the summary honest. +summary["dt"] = parameters["MaxTimeStep"].GetDouble() +summary["final_time"] = parameters["FinalTime"].GetDouble() +with open("results_summary.json", "w") as _f: + json.dump(summary, _f, indent=2) +for _k, _v in sorted(summary.items()): + print("%s=%s" % (_k, _v)) + +# Kratos DEM deletes particles that leave the bounding box and says nothing -- +# an unstable time step therefore looks like a clean, successful run. Refuse to +# exit 0 on a deck that quietly lost its model. +if summary["particles_at_end"] != n_expected: + raise SystemExit( + "DEM lost particles: %d of %d remain. The usual cause is MaxTimeStep " + "above the stability limit; lower dt_safety_factor." + % (summary["particles_at_end"], n_expected)) +if summary["wall_conditions"] == 0: + raise SystemExit( + "no wall conditions were read -- check that " + "DEM_FEM_boundary.mdpa exists next to this script; " + "a missing DEM mdpa is only an INFO line, never an error.") +print("DEM analysis completed.") """ +INPUT_PY = (INPUT_PY.replace("__THREADS__", str(omp_threads)) + .replace("__N_PARTICLES__", repr(n_particles))) with open("input.py", "w") as f: - f.write(input_py) + f.write(INPUT_PY) -print(f"Generated Kratos DEM input files:") -print(f" particles.mdpa: {{n_particles}} particles") -print(f" ProjectParametersDEM.json: dt={{dt}}, T={{T_end}}") -print(f" MaterialsDEM.json: Hertz-viscous-Coulomb") -print(f" input.py: DEMAnalysisStage entry point") +print("Kratos DEM deck written:") +print(" %sDEM.mdpa %d CylinderParticle2D particles" + % (problem_name, n_particles)) +print(" %sDEM_FEM_boundary.mdpa %d RigidFace3D3N conditions (%s)" + % (problem_name, cond_id, ", ".join(walls))) +print(" MaterialsDEM.json materials + material_relations + table") +print(" ProjectParametersDEM.json dt=%.3e (%.3g x Rayleigh %.3e), T=%g, %d steps" + % (dt, dt / dt_rayleigh, dt_rayleigh, T_end, n_steps)) +print(" input.py run it with: python input.py") ''' + return header + body KNOWLEDGE = { @@ -237,24 +401,117 @@ def _dem_2d_kratos(params: dict) -> str: "bonding": ["DEM_KDEM (parallel bond)", "DEM_parallel_bond", "DEM_Dempack"], }, "solver_types": ["explicit (velocity Verlet, default)", "explicit (symplectic Euler)"], + + # EVIDENCE BASIS for the 'dem' entries below. Every pitfall marked + # "(Verified by execution 2026-08-07)" was produced by running Kratos + # 10.4.3 (Release, GCC-15.2, OpenMP) in an environment carrying all 28 + # applications, driving DEMAnalysisStage on a shipped DEMApplication + # fixture and injecting exactly one fault per run. Entries marked + # "(Verified from Kratos source ...)" were read off the C++/Python on + # disk and NOT executed; they are labelled individually. + # + # ProjectParametersDEM.json is a FLAT schema and shares almost nothing + # with the FEM ProjectParameters.json used by StructuralMechanics / + # FluidDynamics. Working minimum, all at top level: "problem_name", + # "FinalTime", "MaxTimeStep", "Dimension", "GravityX/Y/Z", + # "ElementType", "TranslationalIntegrationScheme", + # "RotationalIntegrationScheme", "BoundingBox*", plus + # "solver_settings": {"strategy": ..., "material_import_settings": ...}. + # There is NO "problem_data" block and NO "time_stepping" block. + "project_parameters_schema": { + "shape": "flat, top-level keys \u2014 NOT the FEM problem_data/solver_settings layout", + "problem_name": "top level; the four mdpa filenames are built from it", + "FinalTime": "top level; end time. Overrides problem_data.end_time if both exist", + "MaxTimeStep": "top level; the ONLY key that sets dt. Used verbatim", + "solver_settings.strategy": ( + "python module basename inside DEMApplication: 'sphere_strategy', " + "'continuum_sphere_strategy', 'verlet_continuum_sphere_strategy', " + "'ice_continuum_sphere_strategy'. NOT a solver_type string" + ), + "solver_settings.material_import_settings.materials_filename": "MaterialsDEM.json", + }, + "mdpa_files": ( + "DEM reads FOUR separate mdpa files whose names are problem_name " + "concatenated with a fixed tag and no separator: DEM.mdpa " + "(spheres), DEM_FEM_boundary.mdpa (walls), " + "DEM_Clusters.mdpa, DEM_Inlet.mdpa. " + "solver_settings.model_import_settings.input_filename does NOT name them " + "and has no effect on which files are read." + ), + "materials_schema": ( + "MaterialsDEM.json has three mandatory top-level keys and is NOT the FEM " + "'properties' schema: 'materials' (list of {material_name, material_id, " + "Variables}), 'material_relations' (list of {material_names_list, " + "material_ids_list, Variables}) and 'material_assignation_table' (list of " + "[submodelpart_name, material_name_or_id]). Per-particle data " + "(PARTICLE_DENSITY, YOUNG_MODULUS, POISSON_RATIO) goes in 'materials'; " + "per-CONTACT-PAIR data (STATIC_FRICTION, DYNAMIC_FRICTION, " + "COEFFICIENT_OF_RESTITUTION, DEM_DISCONTINUUM_CONSTITUTIVE_LAW_NAME, " + "DEM_CONTINUUM_CONSTITUTIVE_LAW_NAME, DEM_ROLLING_FRICTION_MODEL_NAME) goes " + "in 'material_relations'. There is no 'constitutive_law': {'name': ...} key." + ), "pitfalls": [ - '[Numerical] Time step must satisfy CFL: dt < min(2*sqrt(m/k)) — use DEM_timestep_safety_factor: 0.5 ' - "Signal: solver reports 'Convergence is not achieved' / 'iteration count exceeded' / oscillating residual; reported quantity disagrees with analytic reference by an order-of-magnitude factor.", - '[Numerical] Bounding box must enclose all particles at all times (AutomaticBoundingBoxOption or manual) ' - "Signal: solver reports 'Convergence is not achieved' / 'iteration count exceeded' / oscillating residual; reported quantity disagrees with analytic reference by an order-of-magnitude factor.", - "[Integration] MaterialsDEM.json constitutive_law name must match exactly (e.g. 'DEM_D_Hertz_viscous_Coulomb') " - "Signal: RuntimeError 'KeyError' from JSON parsing OR 'SubModelPart not found' / 'Property ID ... missing' during AnalysisStage.Initialize.", - '[Integration] MDPA uses SphericParticle3D even for 2D problems (Kratos DEM is always 3D internally) ' - "Signal: RuntimeError 'KeyError' from JSON parsing OR 'SubModelPart not found' / 'Property ID ... missing' during AnalysisStage.Initialize.", - '[Integration] Wall geometry: use RigidFace3D3N elements or the walls_process_list in ProjectParameters ' - "Signal: RuntimeError 'KeyError' from JSON parsing OR 'SubModelPart not found' / 'Property ID ... missing' during AnalysisStage.Initialize.", - '[Numerical] For >10k particles: enable MPI via parallel_type: MPI ' - "Signal: solver reports 'Convergence is not achieved' / 'iteration count exceeded' / oscillating residual; reported quantity disagrees with analytic reference by an order-of-magnitude factor.", - '[Integration] RADIUS must be set per-particle in the MDPA or via a process ' - "Signal: RuntimeError 'KeyError' from JSON parsing OR 'SubModelPart not found' / 'Property ID ... missing' during AnalysisStage.Initialize.", - '[Numerical] Output: VTK with sphere glyphs, or GiD .post.bin format ' - "Signal: solver reports 'Convergence is not achieved' / 'iteration count exceeded' / oscillating residual; reported quantity disagrees with analytic reference by an order-of-magnitude factor.", - ], + "[Input] ProjectParametersDEM.json is a FLAT schema; writing the FEM layout (a 'problem_data' block plus 'solver_settings': {'solver_type': ...}) does not work. DEM looks for solver_settings.strategy, which names a Python module inside DEMApplication ('sphere_strategy'), not a solver_type label. Signal: RuntimeError 'Error: Getting a value that does not exist. entry string : strategy' raised from kratos/sources/kratos_parameters.cpp:426 via DEM_analysis_stage.SetSolverStrategy, before any mesh is read. (Verified by execution 2026-08-07.)", + "[Input] DEM builds four mdpa filenames by concatenating problem_name with the fixed tags DEM, DEM_FEM_boundary, DEM_Clusters and DEM_Inlet \u2014 so problem_name 'mycase' means the spheres live in 'mycaseDEM.mdpa', not 'mycase.mdpa'. A missing file is NOT an error: the reader prints one INFO line and returns, and the run completes normally on an empty model. Signal: 'DEM: Input file DEM.mdpa not found. Continuing.' on stdout, then a full successful run whose SpheresPart ends with 0 elements and 0 nodes; measured 0/0 against 4/4 for the same case with the file present. (Verified by execution 2026-08-07.)", + "[Input] The wall file is silently optional in exactly the same way, which is the more dangerous half: particles then fall through the boundary that the user believes exists. Signal: 'DEM: Input file DEM_FEM_boundary.mdpa not found. Continuing.' and RigidFacePart holds 0 conditions while SpheresPart still holds its 4 particles \u2014 measured 0 vs 18 conditions against the same case with the wall file present. No exception, exit code 0. (Verified by execution 2026-08-07.)", + "[Input] solver_settings.model_import_settings.input_filename is inert for DEM. Setting it to the mdpa base name changes nothing, because file paths come from problem_name; its only reader is the restart utility. Signal: a run whose input_filename points at a real, correctly named mdpa still reports 'Input file DEM.mdpa not found. Continuing.' whenever problem_name disagrees with the file on disk. (Verified from Kratos source 10.4.3, DEM_analysis_stage.py GetInputFilePath/GetProblemNameWithPath \u2014 the inert-key half was not separately executed.)", + "[Input] MaterialsDEM.json is NOT the FEM materials schema. A file shaped {'properties': [{'model_part_name':..., 'Material': {'Variables':..., 'constitutive_law': {'name':...}}}]} \u2014 the StructuralMechanics layout \u2014 has none of the keys DEM reads. Signal: RuntimeError 'Error: Getting a value that does not exist. entry string : materials' from materials_assignation_utility.py; dropping only the assignation table instead gives the same error with 'entry string : material_assignation_table'. (Verified by execution 2026-08-07.)", + "[Input] PARTICLE_FRICTION is not a Kratos variable at all \u2014 it appears in zero files of the installed distribution, including the compiled libraries. Friction is a per-CONTACT-PAIR property named STATIC_FRICTION and DYNAMIC_FRICTION, set in material_relations, not in materials. Signal: RuntimeError 'Error: Value type for \"PARTICLE_FRICTION\" not defined' from read_materials_utility while reading MaterialsDEM.json; the same name passed to KratosGlobals.GetVariable raises ValueError 'Kernel.GetVariable() ERROR: Variable PARTICLE_FRICTION is unknown.'. (Verified by execution 2026-08-07.)", + "[Numerical] Particle density is PARTICLE_DENSITY. DENSITY is also a valid Kratos variable, so writing it is accepted by the materials reader and then never read \u2014 Properties[PARTICLE_DENSITY] silently returns the inserted default 0.0, particle mass becomes 0, the integrator divides by it, and the resulting non-finite coordinates fail the bounding-box test so every particle is erased. Signal: the run exits 0 and reports ANALYSIS COMPLETED while SpheresPart goes from 4 elements to 0 elements / 0 nodes; probing Properties shows PARTICLE_DENSITY = 0.0 alongside DENSITY = 4000.0. No warning of any kind. (Verified by execution 2026-08-07.)", + "[Numerical] Omitting YOUNG_MODULUS behaves the same way: the property read inserts 0.0 rather than raising, giving zero contact stiffness and mutually transparent particles. Signal: Properties[YOUNG_MODULUS] reads back exactly 0.0 after a clean Initialize() and the run proceeds to completion with no message. (Verified by execution 2026-08-07.)", + "[Numerical] MaxTimeStep is the only key that sets dt and it is used verbatim \u2014 DEM performs no stability check and prints no critical-time-step warning at any value. The keys that look like safety nets are dead: DeltaTimeSafetyFactor is stored and never read, AutomaticTimestep is never read at all, and DEM_timestep_safety_factor is not a Kratos key (zero occurrences in the distribution). Signal: MaxTimeStep = 1.0 s on a case whose stable step is 5e-4 prints 'DEM: Total number of time steps expected in the calculation: 0', sets DELTA_TIME to 1.0 and jumps straight to the end time \u2014 no error, no warning. (Verified by execution 2026-08-07.)", + "[Input] FinalTime, not problem_data.end_time, decides when a DEM run stops, and the DEM default silently overwrites whatever end_time the user wrote. Signal: a deck with 'problem_data': {'end_time': 10.0} and no FinalTime key prints 'DEM: Total number of time steps expected in the calculation: 100' and stops at TIME = 0.05 \u2014 the FinalTime default, 200x short \u2014 with no warning. (Verified by execution 2026-08-07.)", + "[Input] Every sub model part that carries entities needs a row in material_assignation_table, including the wall part. Omitting the wall row while the wall mdpa is present is a hard error raised deep in the solver, not at materials-reading time. Signal: RuntimeError 'PROPERTIES_ID is not set for SubModelPart' '. Make sure the Materials file contains material assignation for this SubModelPart', from explicit_solver_strategy InitializeFEMElements and wrapped at runtime with an 'Error: ' prefix. The two quoted halves are separate literals in libKratosDEMCore.so with an interpolation slot between them, so what appears in the middle is whatever identifies the offending part on that path: measured here as the sub model part NAME (DEM-FEM-Wall_floor, for a wall read from an mdpa) and reported previously as a numeric id (1). Those are runtime values, not literals, which is why neither is quoted here. Both are real renderings of one message — do not treat either as a correction of the other, and do not grep the concatenated sentence, which matches nothing. (Verified by execution 2026-08-07; fixture dem_every_entity_submodelpart_needs_an_assignation_row.)", + "[Input] Unknown or misspelled top-level keys are rejected by ValidateAndAssignDefaults, and so are unknown keys inside solver_settings \u2014 but validation is flat, so a typo nested inside any other sub-block is silently ignored. Signal: 'MaxTimestep' (lowercase s) raises RuntimeError 'Error: The item with name \"MaxTimestep\" is present in this Parameters but NOT in the default values', followed by a several-hundred-line dump of both the supplied and the default parameter trees; the same happens for an unknown key under solver_settings. (Verified by execution 2026-08-07.)", + "[API] The contact law name is set per contact PAIR, as DEM_DISCONTINUUM_CONSTITUTIVE_LAW_NAME inside a material_relations entry \u2014 not as a 'constitutive_law' block on a material. The lookup is an unchecked dictionary get, so a wrong or abbreviated name never names itself in the error. Signal: an abbreviated name such as 'DEM_D_Hertz' produces TypeError \"'NoneType' object is not callable\" raised inside sphere_strategy.ModifySubProperties; the offending string is never echoed. Registered names include DEM_D_Hertz_viscous_Coulomb, DEM_D_Linear_viscous_Coulomb, DEM_D_Hertz_viscous_Coulomb_JKR/_DMT, DEM_KDEM, DEM_parallel_bond, DEM_Dempack. (Verified by execution 2026-08-07.)", + "[API] Kratos DEM is NOT 3D-only, contrary to a long-standing catalog claim. There is no SphericParticle2D, but genuine 2D particle elements exist under the Cylinder* stem and construct successfully: CylinderParticle2D and CylinderContinuumParticle2D. The right fix for a planar problem is CylinderParticle2D, not a 3D sphere with constrained out-of-plane DOFs. Signal: CreateNewElement('CylinderParticle2D', ...) succeeds on the same model part where CreateNewElement('SphericParticle2D', ...) raises 'The Element \"SphericParticle2D\" is not registered!'; the 2D discontinuum law is spelled DEM_D_Hertz_viscous_Coulomb2D. (Verified by execution 2026-08-07 \u2014 corrects the earlier 'DEM is always 3D internally' wording.)", + "[API] Wall geometry is built from CONDITIONS, not elements, and only 3D face spellings are registered. Signal: CreateNewCondition('RigidFace3D3N' / 'RigidFace3D4N', ...) constructs, while RigidFace2D2N and RigidEdge3D2N raise 'The Condition \"...\" is not registered!' and list the available DEM conditions. (Verified by execution 2026-08-07.)", + "[API] RADIUS is a core KratosMultiphysics nodal variable, not a DEMApplication one, and it is the single mandatory per-node input in the spheres mdpa. Signal: KratosMultiphysics.RADIUS resolves while DEMApplication.RADIUS raises AttributeError; a particle whose RADIUS was never written keeps the silent default of zero rather than raising. (Verified by execution 2026-08-07.)", + "[Output] There is no working VTK path for a DEM SpheresPart on a stock build, and it fails three different ways. (1) DEM particle geometries have no VTK writer at all. Signal: KM.VtkOutput(spheres_part, settings).PrintOutput() raises RuntimeError 'Modelpart contains elements or conditions with geometries for which no VTK-output is implemented!', wrapped at runtime with an Error: prefix and a trailing Cell type: 30, and located at kratos/input_output/vtk_output.cpp:451, VtkOutput::WriteCell — measured for both CylinderParticle2D and SphericParticle3D, which report the same cell type. The prefix and the cell-type suffix are deliberately left OUTSIDE the quotes: they are assembled around the message at runtime, so a quoted string carrying them matches nothing on disk. (2) The core vtk_output_process cannot even get that far, because its Check() runs before DEM has populated the part's ProcessInfo: output_control_type 'step' raises 'Error: STEP not found in process info of SpheresPart.' and 'time' raises the same for TIME, both from OutputController::Check() at kratos/controllers/output_controller.cpp:81. That one is ordering, not absence — DEM_analysis_stage.Initialize does its own work and only then calls AnalysisStage.Initialize; probing the same ProcessInfo after Initialize returns shows TIME present (STEP still absent). (3) post_vtk_option: true raises ModuleNotFoundError('pyevtk') from SetGraphicalOutput, because DEMApplication.dem_vtk_output opens with 'from pyevtk import hl' and pyevtk is neither a Kratos dependency nor shipped with it. What does work is post_gid_option: true, which writes _Post_Files. (Verified by execution 2026-08-07.)", + # GREP TRAP, and the reason this entry was briefly deleted as a + # suspected fabrication on 2026-08-07 before being restored the same + # day. The geometry message above does not exist as one literal + # anywhere in the distribution. libKratosCore.so holds it as TWO + # adjacent literals with a runtime slot between them: + # + # b'Modelpart contains elements or conditions with ' + # b'geometries for which no VTK-output is implemented!' + # + # so a fixed-string search for the WHOLE SENTENCE returns 0 files + # whether you write one space or two between "with" and + # "geometries". It is not invented: the execution above prints it. + # On this build (Kratos 10.4.3 at /mnt/kratos-tier2/kv) the slot + # renders EMPTY, giving one space; a build where it renders + # non-empty shows an extra token there, which is why the same + # message has been reported with two. + # + # BE PRECISE ABOUT WHAT FAILED HERE, because the honest version is + # less flattering than the convenient one. The project's own + # auditor was NOT fooled by this string. audit_quoted_diagnostics + # falls back to longest_found_prefix when a whole fragment misses, + # and the leading run "Modelpart contains elements or conditions + # with" IS a literal, so even the older left-anchored version + # resolved it and never flagged the entry. What produced the wrong + # deletion was a hand-run `grep -r -a -F` of the full sentence, + # returning 0 and being read as proof — going around the guard that + # exists for exactly this. Grep the interior, or use the auditor. + # + # A DIFFERENT and genuinely harder shape sits next to it: every + # KRATOS_ERROR prepends "Error: " at runtime, so for those the + # literal in the binary is an INTERIOR run and NO leading-prefix + # search can reach it -- "Error: STEP not found in process info of + # SpheresPart." scores 0 at every leading slice while + # "not found in process info of" is present. Those need a + # contiguous-window search. PROPERTIES_ID (see the + # material_assignation_table entry) is split the same way. + ], + "guidance": [ + "[Numerical] There is no automatic time step. Size dt by hand from the Rayleigh/Hertz estimate dt ~ 0.34*sqrt(m/(E*pi*R)) and set it in MaxTimeStep; the formula exists in the source (SphericContinuumParticle::Calculate for DELTA_TIME) but is never invoked by any solver path.", + "[Numerical] The one opt-in automatic-dt path is the automatic_dt_process, added by hand to 'processes'; it prints 'Automatic DT process: Calculated critical time step: seconds.'. It casts unconditionally to SphericContinuumParticle, so pointing it at a plain SphericParticle3D model part dereferences a null pointer. (Source-read, not executed.)", + "[Numerical] BoundingBoxOption defaults to false. With it true, particles leaving the box are deleted with no message and no counter, and mass is silently not conserved; the check runs only on neighbour-search steps. With it false nothing is ever deleted and escapees integrate forever.", + "[Numerical] NeighbourSearchFrequency defaults to 50 and SearchTolerance to 0.0, so contacts are rediscovered only every 50 steps with a search radius equal to the particle radius \u2014 fast particles tunnel through each other silently.", + "[Numerical] For >10k particles: enable MPI via parallel_type: MPI.", + ] }, } From cd8011428377c265bfe398ae1209e9922b215b3f Mon Sep 17 00:00:00 2001 From: Alexander Hermann Date: Sat, 8 Aug 2026 01:30:19 +0200 Subject: [PATCH 2/2] Clamp omp_threads to at least 1 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Flagged in review of #53. `omp_threads=0` reached the deck verbatim, and OMP_NUM_THREADS=0 is not "let OpenMP decide" — it is undefined behaviour, with some runtimes aborting and others silently serialising. A user passing 0 to mean "use the default" would get a deck that fails for a reason unrelated to their simulation. Worse in this particular case: a DEM run reports particle loss silently, so a hard crash is the kinder of the two outcomes and a silent one would be the dangerous one. Verified: omp_threads=0 and -1 both emit 1; 2 passes through unchanged. Co-Authored-By: Claude Opus 5 (1M context) --- src/backends/kratos/generators/dem.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/backends/kratos/generators/dem.py b/src/backends/kratos/generators/dem.py index 43be5914..e9304a7c 100644 --- a/src/backends/kratos/generators/dem.py +++ b/src/backends/kratos/generators/dem.py @@ -26,7 +26,13 @@ def _dem_2d_kratos(params: dict) -> str: dt = params.get("dt", None) dt_safety = float(params.get("dt_safety_factor", 0.02)) problem_name = str(params.get("problem_name", "granular")) - threads = int(params.get("omp_threads", 2)) + # Clamp: OMP_NUM_THREADS=0 is not "use the default", it is undefined + # behaviour — some runtimes abort, others silently serialise. A user + # passing omp_threads=0 meaning "let OpenMP decide" would get a deck that + # fails for a reason unrelated to anything they were trying to simulate, + # and the DEM run reports particle loss silently, so a crash here is the + # kinder outcome of the two. Flagged in review of PR #53. + threads = max(1, int(params.get("omp_threads", 2))) header = f'''\ """Kratos DEM — 2D granular settling — writes a complete DEMApplication deck.