add kratos into the same Snakefile as fenics, still WIP due to a single python function

Author: joergfunger

benchmarks/linear-elastic-plate-with-hole/FEniCS/fenics/environment_simulation.yml

@@ -7,7 +7,6 @@ channel_priority: strict
 dependencies:
   - python=3.12
   - fenics-dolfinx=0.9.*
-  - h5py
   - libadios2=2.10.1
   - petsc4py
   - pint

benchmarks/linear-elastic-plate-with-hole/FEniCS/generate_config.py

@@ -46,7 +46,7 @@ def get_configuration(file):
     "result_dir": result_dir,
     "configuration_to_parameter_file": configuration_to_parameter_file,
     "configurations": list(configurations.values()),
-    "tools": ["fenics"],
+    "tools": ["fenics", "kratos"],
     "benchmark": "linear-elastic-plate-with-hole"
 }
 with open("workflow_config.json", "w") as f:

benchmarks/linear-elastic-plate-with-hole/FEniCS/kratos/Snakefile

@@ -0,0 +1,38 @@
+import json
+
+result_dir = config["result_dir"]
+configuration_to_parameter_file = config["configuration_to_parameter_file"]
+configurations = config["configurations"]
+tool = "kratos"
+
+kratos_input_template = "{tool}/input_template.json"
+kratos_material_template = "{tool}/StructuralMaterials_template.json"
+
+rule run_simulation:
+    input: 
+        script = "{tool}/run_simulation.py",
+        parameters = lambda wildcards: configuration_to_parameter_file[wildcards.configuration],
+        mesh = f"{result_dir}/mesh/mesh_{{configuration}}.msh",
+        kratos_input_template = kratos_input_template,
+        kratos_material_template = kratos_material_template,
+    output:
+        mdpa = f"{result_dir}/{{tool}}/mesh_{{configuration}}.mdpa",
+        kratos_inputfile = f"{result_dir}/{{tool}}/ProjectParameters_{{configuration}}.json",
+        kratos_materialfile = f"{result_dir}/{{tool}}/MaterialParameters_{{configuration}}.json",
+        zip = f"{result_dir}/{{tool}}/solution_field_data_{{configuration}}.zip",
+        metrics = f"{result_dir}/{{tool}}/solution_metrics_{{configuration}}.json",
+    conda:
+        "environment_simulation.yml",
+    shell:
+        """
+        python3 {input.script} \
+            --input_parameter_file {input.parameters} \
+            --input_mesh_file {input.mesh} \
+            --input_kratos_input_template {input.kratos_input_template} \
+            --input_material_template {input.kratos_material_template} \
+            --output_mdpa_file {output.mdpa} \
+            --output_kratos_inputfile {output.kratos_inputfile} \
+            --output_kratos_materialfile {output.kratos_materialfile} \
+            --output_solution_file_zip {output.zip} \
+            --output_metrics_file {output.metrics}
+        """

benchmarks/linear-elastic-plate-with-hole/FEniCS/kratos/environment_simulation.yml

@@ -3,9 +3,11 @@ channels:
   - conda-forge
 dependencies:
   - python=3.10
+  - meshio
   - python-gmsh
   - pint
-  - meshio
+  - pyvista
   - pip
   - pip:
     - KratosMultiphysics-all
+  - sympy

benchmarks/linear-elastic-plate-with-hole/FEniCS/kratos/run_simulation.py

@@ -133,7 +133,7 @@ def msh_to_mdpa(parameter_file: str, mesh_file: str, mdpa_file: str):
 
 def create_kratos_input(
     parameter_file: str,
-    mesh_file: str,
+    mdpa_file: str,
     kratos_input_template_file: str,
     kratos_material_template_file: str,
     kratos_input_file: str,
@@ -199,8 +199,9 @@ def create_kratos_input(
 
     with open(kratos_input_template_file) as f:
         project_parameters_string = f.read()
+    #remove the file extension .mdpa
     project_parameters_string = project_parameters_string.replace(
-        r"{{MESH_FILE}}", mesh_file
+        r"{{MESH_FILE}}", os.path.splitext(mdpa_file)[0]
     )
     project_parameters_string = project_parameters_string.replace(
         r"{{MATERIAL_FILE}}", kratos_material_file
@@ -282,9 +283,9 @@ def create_kratos_input(
 
     create_kratos_input(
         parameter_file=args.input_parameter_file,
-        mesh_file=args.input_mesh_file,
-        input_template_file=args.input_kratos_input_template,
-        material_template_file=args.input_kratos_material_template,
+        mdpa_file=args.output_mdpa_file,
+        kratos_input_template_file=args.input_kratos_input_template,  # <-- fixed argument name
+        kratos_material_template_file=args.input_material_template,
         kratos_input_file=args.output_kratos_inputfile,
         kratos_material_file=args.output_kratos_materialfile,
     )
@@ -306,7 +307,7 @@ def create_kratos_input(
         parameters = json.load(f)
     config = parameters["configuration"]
     output_dir = Path(os.path.dirname(os.path.abspath(args.output_kratos_inputfile))) / str(config)
-    mesh = pyvista.read("Structure_0_1.vtk")
+    mesh = pyvista.read(str(output_dir / "Structure_0_1.vtk"))
     max_von_mises_stress = float(mesh["VON_MISES_STRESS"].max())
     print("Max Von Mises Stress:", max_von_mises_stress)
     metrics = {
@@ -315,15 +316,9 @@ def create_kratos_input(
     with open(args.output_metrics_file, "w") as f:
         json.dump(metrics, f, indent=4)
 
-    file_patterns = [
-        str(output_dir / "Structure_0_1.vtk"),
-    ]
+    files_to_store = [str(output_dir / "Structure_0_1.vtk")]
 
-    files_to_store = []
-    for pattern in file_patterns:
-        files_to_store.extend(Path().glob(pattern))
-    
     import zipfile
     with zipfile.ZipFile(args.output_solution_file_zip, "w") as zipf:
         for filepath in files_to_store:
-            zipf.write(filepath, arcname=f"result_{config}.vtk") 
+            zipf.write(filepath, arcname=f"result_{config}.vtk")

benchmarks/linear-elastic-plate-with-hole/FEniCS/plateWithHoleSolution.py

@@ -119,13 +119,14 @@ def displacement_symbolic_str(self, X_str: str, Y_str: str):
         ux, uy = self.displacement_symbolic()
         # Define symbols
         x, y = sp.symbols('x y')
-        # Substitute variable names
-        subs = {x: sp.Symbol(X_str), y: sp.Symbol(Y_str)}
-        ux_sub = ux.subs(subs)
-        uy_sub = uy.subs(subs)
-        # Convert to string
-        ux_str = str(ux_sub)
-        uy_str = str(uy_sub)
+        E, nu, a, T = sp.symbols('E nu a T')
+        subs_vars = {x: sp.Symbol(X_str), y: sp.Symbol(Y_str)}
+        subs_params = {E: self.E, nu: self.nu, a: self.radius, T: self.load}
+        ux_sub = ux.subs(subs_vars).subs(subs_params)
+        uy_sub = uy.subs(subs_vars).subs(subs_params)
+        # Convert to string using sympy.sstr for single-line output with **
+        ux_str = sp.sstr(ux_sub)
+        uy_str = sp.sstr(uy_sub)
         return ux_str, uy_str
 
     def stress(self, x: np.ndarray) -> np.ndarray:
@@ -141,5 +142,3 @@ def stress(self, x: np.ndarray) -> np.ndarray:
         sxy = self.sxy_func(arr[0], arr[1])
         syy = self.syy_func(arr[0], arr[1])
         return sxx, sxy, syy
-
-