diff --git a/.github/workflows/run-benchmark.yml b/.github/workflows/run-benchmark.yml index 3f37170..85e71ea 100644 --- a/.github/workflows/run-benchmark.yml +++ b/.github/workflows/run-benchmark.yml @@ -37,126 +37,18 @@ jobs: - name: Update environment run: mamba env update -n model-validation -f environment_benchmarks.yml - - name: generate-config-files + - name: Run Example linear-elastic-plate-with-hole_fenics using Benchmarking package shell: bash -l {0} run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - python generate_config.py - - - name: run_linear-elastic-plate-with-hole-benchmarks_snakemake - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - snakemake --use-conda --force --cores 'all' - snakemake --use-conda --force --cores all \ - --reporter metadata4ing \ - --report-metadata4ing-paramscript ../common/parameter_extractor.py \ - --report-metadata4ing-config metadata4ing.config \ - --report-metadata4ing-filename $SNAKEMAKE_RESULT_FILE - - - name: run_linear-elastic-plate-with-hole-benchmarks_nextflow - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - nextflow run main.nf -params-file workflow_config.json -c ../common/nextflow.config -plugins nf-prov@1.4.0 + cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/fenics/ + python run_benchmark.py - name: Archive Linear Elastic plate with a hole benchmark data for snakemake uses: actions/upload-artifact@v4 with: name: snakemake_results_linear-elastic-plate-with-hole path: | - benchmarks/linear-elastic-plate-with-hole/${{ env.SNAKEMAKE_RESULT_FILE }}.zip - - - name: Archive Linear Elastic plate with a hole benchmark data for nextflow - uses: actions/upload-artifact@v4 - with: - name: nextflow_results_linear-elastic-plate-with-hole - path: | - benchmarks/linear-elastic-plate-with-hole/nextflow_results/ + examples/linear-elastic-plate-with-hole/fenics/results/ - process-artifacts: - runs-on: ubuntu-latest - needs: run-simulation - steps: - - name: Checkout repo content - uses: actions/checkout@v2 - - - name: Download artifact - uses: actions/download-artifact@v4 - with: - name: snakemake_results_linear-elastic-plate-with-hole - path: ./artifact_files - - name: Unzip Snakemake Result File - run: | - mkdir -p ./$SNAKEMAKE_RESULT_FILE - unzip -o ./artifact_files/$SNAKEMAKE_RESULT_FILE.zip -d ./$SNAKEMAKE_RESULT_FILE - - - name: Setup Mambaforge with postprocessing env - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-version: latest - activate-environment: postprocessing - use-mamba: true - environment-file: benchmarks/linear-elastic-plate-with-hole/environment_postprocessing.yml - - - name: Validate Snakemake Result File - shell: bash -l {0} - run: | - python benchmarks/common/validate_provenance.py \ - --provenance_folderpath "./$SNAKEMAKE_RESULT_FILE" - - - name: Run plotting script - shell: bash -l {0} - run: | - python benchmarks/linear-elastic-plate-with-hole/plot_metrics.py \ - --provenance_folderpath "./$SNAKEMAKE_RESULT_FILE" \ - --output_file $PROVENANACE_FILE_NAME - - - name: Upload snakemake results file as artifact - uses: actions/upload-artifact@v4 - with: - name: element-size-vs-stress-plot - path: ${{ env.PROVENANACE_FILE_NAME }} - - - name: Re-zip snakemake result folder - run: | - cd "./${SNAKEMAKE_RESULT_FILE}" - zip -r "../${SNAKEMAKE_RESULT_FILE}.zip" . - - - name: Upload RoCrate Zip file onto RoHub - id: rohub_upload - shell: bash -l {0} - continue-on-error: true - run: | - LOG_FILE=rohub_upload_error.log - - # Run Python and capture all stdout+stderr - python benchmarks/common/upload_provenance.py \ - --provenance_folderpath "./${SNAKEMAKE_RESULT_FILE}.zip" \ - --benchmark_name "linear-elastic-plate-with-hole" \ - --username "${{ secrets.ROHUB_USERNAME }}" \ - --password "${{ secrets.ROHUB_PASSWORD }}" \ - 2>&1 | tee "$LOG_FILE" - - PYTHON_EXIT_CODE=${PIPESTATUS[0]} - - # Check exit code: 0 = failure, 1 = success (per your convention) - if [ $PYTHON_EXIT_CODE -eq 0 ]; then - echo "=== RoHub upload failed — see log below ===" - cat "$LOG_FILE" - else - echo "✅ RoHub upload succeeded" - rm -f "$LOG_FILE" - fi - - # Export exit code for subsequent steps if needed - echo "python_exit_code=$PYTHON_EXIT_CODE" >> $GITHUB_ENV - - - name: Upload RoHub error log - if: ${{ env.python_exit_code == '0' }} - uses: actions/upload-artifact@v4 - with: - name: rohub-upload-error-log - path: rohub_upload_error.log diff --git a/benchmarks/linear-elastic-plate-with-hole/fenics/fenics.nf b/benchmarks/linear-elastic-plate-with-hole/fenics/fenics.nf deleted file mode 100644 index a655fe8..0000000 --- a/benchmarks/linear-elastic-plate-with-hole/fenics/fenics.nf +++ /dev/null @@ -1,35 +0,0 @@ -params.tool = "fenics" - -process run_simulation { - publishDir "${params.result_dir}/${params.tool}/" - conda './fenics/environment_simulation.yml' - - input: - path python_script - tuple val(configuration), path(parameter_file), path(mesh_file) - - - output: - tuple val(configuration), path("solution_field_data_${configuration}.zip"), path("solution_metrics_${configuration}.json") - - script: - """ - python3 $python_script --input_parameter_file $parameter_file --input_mesh_file $mesh_file --output_solution_file_zip "solution_field_data_${configuration}.zip" --output_metrics_file "solution_metrics_${configuration}.json" - """ -} - -workflow fenics_workflow { - - take: - mesh_data // tuple(configuration, parameters, mesh) - result_dir - - main: - params.result_dir = result_dir - run_sim_script = Channel.value(file('fenics/run_fenics_simulation.py')) - output_process_run_simulation = run_simulation( run_sim_script, mesh_data ) - - emit: - output_process_run_simulation - -} \ No newline at end of file diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/kratos.nf b/benchmarks/linear-elastic-plate-with-hole/kratos/kratos.nf deleted file mode 100644 index c5d8b91..0000000 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/kratos.nf +++ /dev/null @@ -1,123 +0,0 @@ -params.tool = "kratos" - -process mesh_to_mdpa { - publishDir "${params.result_dir}/${params.tool}/" - conda './kratos/environment_simulation.yml' - - input: - path python_script - tuple val(configuration), path(parameter_file), path(mesh_file) - - output: - tuple val(configuration), path("mesh_${configuration}.mdpa") - - script: - """ - python3 ${python_script} \ - --input_parameter_file ${parameter_file} \ - --input_mesh_file ${mesh_file} \ - --output_mdpa_file mesh_${configuration}.mdpa - """ -} - -process create_kratos_input_and_run_simulation { - - // The process combines the creation of the Kratos input file (json) and the execution of the simulation. Initially, these were two separate processes. - // The combination was necessary because the create_kratos_input.py specifies the location of the mesh file and the output location of the simulation results in - // the json file. In the case of Nextflow, these locations are related to the process's sub-directory (inside the work directory). Executing the simulation as a - // separate process results in a failure to find the mesh file and write the output files unless the paths (in the json file) are explicitly provided as an input - // to the process. - // This is not an issue in the case of Snakemake, as the working directory doesn't automatically change between different rules. - - publishDir "${params.result_dir}/${params.tool}/" - conda './kratos/environment_simulation.yml' - - input: - path script_create_kratos_input - path script_run_kratos - tuple val(configuration), path(parameters), path(mdpa) - path kratos_input_template - path kratos_material_template - - output: - tuple val(configuration), path("ProjectParameters_${configuration}.json"), path("MaterialParameters_${configuration}.json"), path("${configuration}/Structure_0_1.vtk") - - script: - """ - python3 ${script_create_kratos_input} \ - --input_parameter_file ${parameters} \ - --input_mdpa_file ${mdpa} \ - --input_kratos_input_template ${kratos_input_template} \ - --input_material_template ${kratos_material_template} \ - --output_kratos_inputfile ProjectParameters_${configuration}.json \ - --output_kratos_materialfile MaterialParameters_${configuration}.json - - python3 ${script_run_kratos} \ - --input_parameter_file ${parameters} \ - --input_kratos_inputfile "ProjectParameters_${configuration}.json" \ - --input_kratos_materialfile "MaterialParameters_${configuration}.json" - """ -} - -process postprocess_kratos_results { - publishDir "${params.result_dir}/${params.tool}/" - conda './kratos/environment_simulation.yml' - - input: - path python_script - tuple val(configuration), path(parameter_file), path(result_vtk) - - output: - tuple val(configuration), path("solution_field_data_${configuration}.zip"), path("solution_metrics_${configuration}.json") - - script: - """ - python3 ${python_script} \ - --input_parameter_file ${parameter_file} \ - --input_result_vtk ${result_vtk} \ - --output_solution_file_zip solution_field_data_${configuration}.zip \ - --output_metrics_file solution_metrics_${configuration}.json - """ -} - -workflow kratos_workflow { - take: - mesh_data // tuple(configuration, parameters, mesh) //change the name - result_dir - - main: - params.result_dir = result_dir - - // Define script paths - msh_to_mdpa_script = Channel.value(file('kratos/msh_to_mdpa.py')) - create_input_script = Channel.value(file('kratos/create_kratos_input.py')) - run_sim_script = Channel.value(file('kratos/run_kratos_simulation.py')) - postprocess_script = Channel.value(file('kratos/postprocess_results.py')) - - // Template files - kratos_input_template = Channel.value(file('kratos/input_template.json')) - kratos_material_template = Channel.value(file('kratos/StructuralMaterials_template.json')) - - // Process pipeline - output_process_mesh_to_mdpa = mesh_to_mdpa(msh_to_mdpa_script, mesh_data) - - input_process_create_kratos_input = mesh_data.join(output_process_mesh_to_mdpa).map { tuple(it[0], it[1], it[3]) } - - //input_process_create_kratos_input.view() - output_create_kratos_input_and_run_simulation = create_kratos_input_and_run_simulation( - create_input_script, - run_sim_script, - input_process_create_kratos_input, - kratos_input_template, - kratos_material_template - ) - - input_process_postprocess_kratos_results = mesh_data.join(output_create_kratos_input_and_run_simulation).map { tuple(it[0], it[1], it[5]) } - - - output_process_postprocess_kratos_results = postprocess_kratos_results(postprocess_script,input_process_postprocess_kratos_results) - - emit: - output_process_postprocess_kratos_results -} - diff --git a/benchmarks/linear-elastic-plate-with-hole/main.nf b/benchmarks/linear-elastic-plate-with-hole/main.nf deleted file mode 100644 index e6bf069..0000000 --- a/benchmarks/linear-elastic-plate-with-hole/main.nf +++ /dev/null @@ -1,173 +0,0 @@ - -include { fenics_workflow } from './fenics/fenics.nf' -include { kratos_workflow } from './kratos/kratos.nf' - -process create_mesh { - //publishDir "$result_dir/mesh/" - publishDir "${params.result_dir}/mesh/" - conda 'environment_mesh.yml' - - input: - path python_script - val configuration - path parameter_file - - output: - // val(configuration) works as matching key with the input channel in the workflow - tuple val(configuration), path("mesh_${configuration}.msh") - - script: - """ - python3 $python_script --input_parameter_file $parameter_file --output_mesh_file "mesh_${configuration}.msh" - """ -} - -process summary{ - publishDir "${params.result_dir}/${tool}/" - conda 'environment_postprocessing.yml' - - input: - path python_script - val configuration - val parameter_file - val mesh_file - val solution_metrics - val solution_field_data - val benchmark - val benchmark_uri - val tool - - output: - path("summary.json") - - script: - """ - python3 $python_script \ - --input_configuration ${configuration.join(' ')} \ - --input_parameter_file ${parameter_file.join(' ')} \ - --input_mesh_file ${mesh_file.join(' ')} \ - --input_solution_metrics ${solution_metrics.join(' ')} \ - --input_solution_field_data ${solution_field_data.join(' ')} \ - --input_benchmark ${benchmark} \ - --input_benchmark_uri ${benchmark_uri} \ - --output_summary_json "summary.json" - - """ -} - - -def prepare_inputs_for_process_summary(input_process_run_simulation, output_process_run_simulation) { - - // Input: channels of the input and the output of the simulation process - // Output: a tuple of channels to be used as input for the summary process - // Purpose: To prepare inputs for the summary process (invoked once per simulation tool) from the output of the simulation process (depending on the number of configurations, invoked multiple times per simulation tool). - - // Firstly, the join operation is performed between the input and output channels of the simulation process based on a matching key (configuration). - - // Secondly, the individual components (configuration, parameter_file, mesh_file, solution_field_data, solution_metrics) are extracted from the joined tuples and collected into separate lists. - // The collect() method outputs a channel with a single entry as the summary process runs once per simulation tool. This single entry is a list of all the configurations or parameter files or mesh files etc. - - def matched_channels = input_process_run_simulation.join(output_process_run_simulation) - - def branched_channels = matched_channels.multiMap{ v, w, x, y, z -> - configuration : v - parameter_file : w - mesh : x - solution_field : y - metrics : z } - - return [ - branched_channels.configuration.collect(), - branched_channels.parameter_file.collect(), - branched_channels.mesh.collect(), - branched_channels.solution_field.collect(), - branched_channels.metrics.collect() - ] -} - -workflow { - main: - - def parameter_files_path = [] - params.configurations.each { elem -> - parameter_files_path.add(file(params.configuration_to_parameter_file[elem])) - } - - def ch_parameter_files = Channel.fromList(parameter_files_path) - def ch_configurations = Channel.fromList(params.configurations) - def ch_mesh_python_script = Channel.value(file('create_mesh.py')) - - //Creating Mesh - - output_process_create_mesh = create_mesh(ch_mesh_python_script, ch_configurations, ch_parameter_files) - - input_process_run_simulation = ch_configurations.merge(ch_parameter_files).join(output_process_create_mesh) - - //Running Simulation - - ch_tools = Channel.fromList(params.tools) - - input_process_run_simulation_with_tool = ch_tools.combine(input_process_run_simulation) - input_fenics_workflow = input_process_run_simulation_with_tool.filter{ it[0] == 'fenics' }.map{_w,x,y,z -> tuple(x,y,z)} - input_kratos_workflow = input_process_run_simulation_with_tool.filter{ it[0] == 'kratos' }.map{_w,x,y,z -> tuple(x,y,z)} - - - fenics_workflow(input_fenics_workflow, params.result_dir) - output_fenics_workflow = fenics_workflow.out - def (fenics_configurations,\ - fenics_parameter_files,\ - fenics_meshes,\ - fenics_solution_fields,\ - fenics_summary_metrics) = prepare_inputs_for_process_summary(input_fenics_workflow, output_fenics_workflow) - - - kratos_workflow(input_kratos_workflow, params.result_dir) - output_kratos_workflow = kratos_workflow.out - def (kratos_configurations, \ - kratos_parameter_files, \ - kratos_meshes, \ - kratos_solution_fields, \ - kratos_summary_metrics) = prepare_inputs_for_process_summary(input_kratos_workflow, output_kratos_workflow) - - - // channels are concatenated in the same order as they are passed to the .concat. The order should be consistent with the order of tools in ch_tools. - input_summary_configuration = fenics_configurations.concat(kratos_configurations) - input_summary_parameter_file = fenics_parameter_files.concat(kratos_parameter_files) - input_summary_mesh = fenics_meshes.concat(kratos_meshes) - input_summary_solution_field = fenics_solution_fields.concat(kratos_solution_fields) - input_summary_metrics = fenics_summary_metrics.concat(kratos_summary_metrics) - - //Summarizing results - def ch_benchmark = Channel.value(params.benchmark) - def ch_benchmark_uri = Channel.value(params.benchmark_uri) - def ch_summarize_python_script = Channel.value(file('../common/summarize_results.py')) - summary(ch_summarize_python_script, \ - input_summary_configuration, \ - input_summary_parameter_file, \ - input_summary_mesh, \ - input_summary_metrics, \ - input_summary_solution_field, \ - ch_benchmark, \ - ch_benchmark_uri, \ - ch_tools) - -} -/* -Steps to add a new simulation tool to the workflow: - -1. Write the tool-specific workflow, scripts, environment file and store them in the benchmarks/linear-elastic-plate-with-hole/tool_name/. -2. Add the tool name to "tools" workflow_config.json (generated here using generate_config.py) -3. Include the tool-specific workflow script at the top of this file. -4. Create an input channel for the new tool (e.g. see the definition of input_fenics_workflow) -5. Invoke the new tool-specific workflow (similar to fenics_workflow) & using its output, prepare inputs for the summary process. -6. Concatenate the prepared inputs to form the final input channels for the summary process. - ---------------------------------------------------------------------------------------------------------------------------------- - -Remark: Care should be taken to track the entries in the I/O channels, as the process output for a given configuration -may not arrive in the same order as the inputs were sent. When reusing channel entries after process execution, outputs should -be matched with their corresponding inputs using a common key. - -Information on channel operations: https://www.nextflow.io/docs/latest/reference/operator.html -Information on channels: https://training.nextflow.io/2.2/basic_training/channels/ -*/ \ No newline at end of file diff --git a/environment_benchmarks.yml b/environment_benchmarks.yml index 0e36b3a..76feea0 100644 --- a/environment_benchmarks.yml +++ b/environment_benchmarks.yml @@ -4,11 +4,11 @@ channels: - bioconda dependencies: - snakemake - - nextflow - pyvista - meshio - conda-lock - conda - pip - pip: - - "git+https://github.com/izus-fokus/snakemake-report-plugin-metadata4ing@v1.2.3#egg=snakemake-report-plugin-metadata4ing" \ No newline at end of file + - "git+https://github.com/izus-fokus/snakemake-report-plugin-metadata4ing@v1.2.5#egg=snakemake-report-plugin-metadata4ing" + - -e . \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/fenics/environment_simulation.yml b/examples/linear-elastic-plate-with-hole/fenics/environment_simulation.yml new file mode 100644 index 0000000..6ed60ef --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/fenics/environment_simulation.yml @@ -0,0 +1,16 @@ +name: fenics_simulation +# Environment file for fenics simulation scripts. Called by fenics tool workflow. + +channels: + - conda-forge + +channel_priority: strict + +dependencies: + - python=3.12 + - fenics-dolfinx=0.9.* + - mpich + - petsc4py + - pint + - python-gmsh + - sympy diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py new file mode 100644 index 0000000..318c001 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py @@ -0,0 +1,63 @@ +from pathlib import Path +import zipfile +import json +from benchmarking import benchmark + +""" +The script performs the following steps: + +1. Extracts the benchmark files from a zip archive (currently assuming that it is an RO-Crate of the benchmark). +2. Creates a benchmark object and registers the simulation tool (FEniCS) along with its version and URI. +3. Adds the simulation script and environment file to the benchmark object. +4. Iterates through the parameter configuration files, checks the "element-size" value, and if it meets the specified condition (>= 0.025), it generates and executes the Snakemake workflow for that configuration. + +The results of each run (and the files used by it) are stored in the directory with the configuration name. +""" + + +#################################################################################################### +#################################################################################################### +# Benchmark Extraction +#################################################################################################### +#################################################################################################### + +zipped_benchmark_dir = Path(__file__).resolve().parent.parent +unzipped_benchmark_dir = Path(__file__).resolve().parent + +with zipfile.ZipFile(zipped_benchmark_dir / "linear-elastic-plate-with-hole.zip", 'r') as zip_ref: + # Extract all files + zip_ref.extractall(unzipped_benchmark_dir) + + +#################################################################################################### +#################################################################################################### +# Creation of benchmark object and the addition of simulation tool scripts +#################################################################################################### +#################################################################################################### + +linear_elastic_problem = benchmark("linear-elastic-plate-with-hole", \ + benchmark_dir=unzipped_benchmark_dir, \ + benchmark_uri="https://portal.mardi4nfdi.de/wiki/Model:6775296") + +linear_elastic_problem.add_tool("fenics", \ + version="0.9.0", \ + uri="https://github.com/FEniCS/dolfinx") + +linear_elastic_problem.add_tool_scripts( + simulation_script = unzipped_benchmark_dir / "run_fenics_simulation.py", + environment_file = unzipped_benchmark_dir / "environment_simulation.yml" +) + +#################################################################################################### +#################################################################################################### +# Conditional creation and execution of snakemake workflows based on parameter configurations +#################################################################################################### +#################################################################################################### + +for file in unzipped_benchmark_dir.glob("parameters_*.json"): + with open(file, "r") as f: + data = json.load(f) + if data.get("element-size").get("value") >= 0.025: + linear_elastic_problem.generate_workflow(file.name, data.get("configuration")) + linear_elastic_problem.run_workflow() + \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py b/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py new file mode 100644 index 0000000..ef980b2 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py @@ -0,0 +1,301 @@ +import json +import sys +from argparse import ArgumentParser + +from pathlib import Path +import dolfinx as df +import basix.ufl +import numpy as np +import ufl +from dolfinx.fem.petsc import LinearProblem +from petsc4py.PETSc import ScalarType +from mpi4py import MPI +from pint import UnitRegistry + +# Add parent directory to sys.path +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from analytical_solution import AnalyticalSolution + + +def run_fenics_simulation( + parameter_file: str, mesh_file: str, solution_file_zip: str, metrics_file: str +) -> None: + ureg = UnitRegistry() + with open(parameter_file) as f: + parameters = json.load(f) + + mesh, cell_tags, facet_tags = df.io.gmshio.read_from_msh( + mesh_file, + comm=MPI.COMM_WORLD, + gdim=2, + ) + + V = df.fem.functionspace(mesh, ("CG", parameters["element-degree"], (2,))) + + tags_left = facet_tags.find(1) + tags_bottom = facet_tags.find(2) + tags_right = facet_tags.find(3) + tags_top = facet_tags.find(4) + + # Boundary conditions + dofs_left = df.fem.locate_dofs_topological(V.sub(0), 1, tags_left) + dofs_bottom = df.fem.locate_dofs_topological(V.sub(1), 1, tags_bottom) + dofs_right = df.fem.locate_dofs_topological(V, 1, tags_right) + dofs_top = df.fem.locate_dofs_topological(V, 1, tags_top) + + bc_left = df.fem.dirichletbc(0.0, dofs_left, V.sub(0)) + bc_bottom = df.fem.dirichletbc(0.0, dofs_bottom, V.sub(1)) + + E = ( + ureg.Quantity( + parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"] + ) + .to_base_units() + .magnitude + ) + nu = ( + ureg.Quantity( + parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"] + ) + .to_base_units() + .magnitude + ) + radius = ( + ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + .to_base_units() + .magnitude + ) + L = ( + ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + .to_base_units() + .magnitude + ) + load = ( + ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"]) + .to_base_units() + .magnitude + ) + + analytical_solution = AnalyticalSolution( + E=E, + nu=nu, + radius=radius, + L=L, + load=load, + ) + + def eps(v): + return ufl.sym(ufl.grad(v)) + + def sigma(v): + # plane stress + epsilon = eps(v) + return ( + E + / (1.0 - nu**2) + * ((1.0 - nu) * epsilon + nu * ufl.tr(epsilon) * ufl.Identity(2)) + ) + + def as_tensor(v): + return ufl.as_matrix([[v[0], v[2]], [v[2], v[1]]]) + + dx = ufl.Measure( + "dx", + metadata={ + "quadrature_degree": parameters["quadrature-degree"], + "quadrature_scheme": parameters["quadrature-rule"], + }, + ) + ds = ufl.Measure( + "ds", + domain=mesh, + subdomain_data=facet_tags, + ) + stress_space = df.fem.functionspace( + mesh, ("CG", parameters["element-degree"], (2, 2)) + ) + stress_function = df.fem.Function(stress_space) + + u = df.fem.Function(V, name="u") + u_prescribed = df.fem.Function(V, name="u_prescribed") + u_prescribed.interpolate(lambda x: analytical_solution.displacement(x)) + u_prescribed.x.scatter_forward() + + u_ = ufl.TestFunction(V) + v_ = ufl.TrialFunction(V) + a = df.fem.form(ufl.inner(sigma(u_), eps(v_)) * dx) + + # set rhs to zero + f = df.fem.form(ufl.inner(df.fem.Constant(mesh, np.array([0.0, 0.0])), u_) * ufl.ds) + + bc_right = df.fem.dirichletbc(u_prescribed, dofs_right) + bc_top = df.fem.dirichletbc(u_prescribed, dofs_top) + solver = LinearProblem( + a, + f, + bcs=[bc_left, bc_bottom, bc_right, bc_top], + u=u, + petsc_options={ + "ksp_type": "gmres", + "ksp_rtol": 1e-14, + "ksp_atol": 1e-14, + }, + ) + solver.solve() + + def project( + v: df.fem.Function | ufl.core.expr.Expr, + V: df.fem.FunctionSpace, + dx: ufl.Measure = ufl.dx, + ) -> df.fem.Function: + """ + Calculates an approximation of `v` on the space `V` + + Args: + v: The expression that we want to evaluate. + V: The function space on which we want to evaluate. + dx: The measure that is used for the integration. This is important, if + either `V` is a quadrature space or `v` is a ufl expression containing a quadrature space. + + Returns: + A function if `u` is None, otherwise `None`. + + """ + dv = ufl.TrialFunction(V) + v_ = ufl.TestFunction(V) + a_proj = ufl.inner(dv, v_) * dx + b_proj = ufl.inner(v, v_) * dx + + solver = LinearProblem(a_proj, b_proj) + uh = solver.solve() + return uh + + plot_space_stress = df.fem.functionspace( + mesh, ("DG", parameters["element-degree"] - 1, (2, 2)) + ) + plot_space_mises = df.fem.functionspace( + mesh, ("DG", parameters["element-degree"] - 1, (1,)) + ) + stress_nodes_red = project(sigma(u), plot_space_stress, dx) + stress_nodes_red.name = "stress" + + def mises_stress(u): + stress = sigma(u) + p = ufl.tr(stress) / 3.0 + s = stress - p * ufl.Identity(2) + return ufl.as_vector([(3.0 / 2.0) ** 0.5 * (ufl.inner(s, s) + p * p) ** 0.5]) + + mises_stress_nodes = project(mises_stress(u), plot_space_mises, dx) + mises_stress_nodes.name = "von_mises_stress" + + # Write each function to its own VTK file on all ranks + output_dir = Path(solution_file_zip).parent + with df.io.VTKFile( + MPI.COMM_WORLD, + str( + output_dir + / f"solution_field_data_displacements_{parameters['configuration']}.vtk" + ), + "w", + ) as vtk: + vtk.write_function(u, 0.0) + with df.io.VTKFile( + MPI.COMM_WORLD, + str( + output_dir / f"solution_field_data_stress_{parameters['configuration']}.vtk" + ), + "w", + ) as vtk: + vtk.write_function(stress_nodes_red, 0.0) + with df.io.VTKFile( + MPI.COMM_WORLD, + str( + output_dir + / f"solution_field_data_mises_stress_{parameters['configuration']}.vtk" + ), + "w", + ) as vtk: + vtk.write_function(mises_stress_nodes, 0.0) + + # extract maximum von Mises stress + max_mises_stress_nodes = np.max(mises_stress_nodes.x.array) + + # Compute von Mises stress at quadrature (Gauss) points and extract maximum (global across MPI) + quad_element = basix.ufl.quadrature_element( + mesh.topology.cell_name(), + value_shape=(1,), + degree=parameters["quadrature-degree"], + ) + + Q_mises = df.fem.functionspace(mesh, quad_element) + mises_qp = df.fem.Function(Q_mises, name="von_mises_stress_qp") + expr_qp = df.fem.Expression(mises_stress(u), Q_mises.element.interpolation_points()) + mises_qp.interpolate(expr_qp) + max_mises_stress_gauss_points = MPI.COMM_WORLD.allreduce( + np.max(mises_qp.x.array), op=MPI.MAX + ) + # Save metrics + metrics = { + "max_von_mises_stress_nodes": max_mises_stress_nodes, + "max_von_mises_stress_gauss_points": max_mises_stress_gauss_points, + } + + if MPI.COMM_WORLD.rank == 0: + with open(metrics_file, "w") as f: + json.dump(metrics, f, indent=4) + # store all .vtu, .pvtu and .vtk files for this configuration in the zip file + import zipfile + + config = parameters["configuration"] + file_patterns = [ + str(output_dir / f"solution_field_data_displacements_{config}*"), + str(output_dir / f"solution_field_data_stress_{config}*"), + str(output_dir / f"solution_field_data_mises_stress_{config}*"), + ] + + files_to_store = [] + for pattern in file_patterns: + files_to_store.extend( + filter( + # filter for all file endings because this is not possible with glob + lambda path: path.suffix in [".vtk", ".vtu", ".pvtu"], + Path().glob(pattern), + ) + ) + # files_to_store.extend(Path().glob(pattern)) + with zipfile.ZipFile(solution_file_zip, "w") as zipf: + for filepath in files_to_store: + zipf.write(filepath, arcname=filepath.name) + + +if __name__ == "__main__": + parser = ArgumentParser( + description="Run FEniCS simulation for a plate with a hole.\n" + "Inputs: --input_parameter_file, --input_mesh_file\n" + "Outputs: --output_solution_file_hdf5, --output_metrics_file" + ) + parser.add_argument( + "--input_parameter_file", + required=True, + help="JSON file containing simulation parameters (input)", + ) + parser.add_argument( + "--input_mesh_file", required=True, help="Path to the mesh file (input)" + ) + parser.add_argument( + "--output_solution_file_zip", + required=True, + help="Path to the zipped solution files (output)", + ) + parser.add_argument( + "--output_metrics_file", + required=True, + help="Path to the output metrics JSON file (output)", + ) + args, _ = parser.parse_known_args() + run_fenics_simulation( + args.input_parameter_file, + args.input_mesh_file, + args.output_solution_file_zip, + args.output_metrics_file, + ) diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/MainKratos.py b/examples/linear-elastic-plate-with-hole/kratos/kratos/MainKratos.py new file mode 100755 index 0000000..02d2394 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/MainKratos.py @@ -0,0 +1,19 @@ +from __future__ import print_function, absolute_import, division #makes KratosMultiphysics backward compatible with python 2.6 and 2.7 + +import KratosMultiphysics +from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis +import sys +""" +For user-scripting it is intended that a new class is derived +from StructuralMechanicsAnalysis to do modifications +""" +parameters = sys.argv[1] + +if __name__ == "__main__": + + with open(parameters,'r') as parameter_file: + parameters = KratosMultiphysics.Parameters(parameter_file.read()) + + model = KratosMultiphysics.Model() + simulation = StructuralMechanicsAnalysis(model,parameters) + simulation.Run() diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/Snakefile b/examples/linear-elastic-plate-with-hole/kratos/kratos/Snakefile new file mode 100644 index 0000000..3e59279 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/Snakefile @@ -0,0 +1,71 @@ +rule mesh_to_mdpa: + input: + parameters = lambda wildcards: configuration_to_parameter_file[wildcards.configuration], + mesh = f"{result_dir}/mesh/mesh_{{configuration}}.msh", + script = f"{tool}/msh_to_mdpa.py", + output: + mdpa = f"{result_dir}/mesh_{{configuration}}.mdpa", + conda: + f"{tool}/environment_simulation.yml", + shell: + """ + python3 {input.script} \ + --input_parameter_file {input.parameters} \ + --input_mesh_file {input.mesh} \ + --output_mdpa_file {output.mdpa} + """ + +rule create_kratos_input_and_run_simulation: + # The process combines the creation of the Kratos input file (json) and the execution of the simulation. Initially, these were two separate processes. + # The combination was necessary because the create_kratos_input.py specifies the location of the mesh file and the output location of the simulation results in + # the json file. In the case of Nextflow, these locations are related to the process's sub-directory (inside the work directory). Executing the simulation as a + # separate process results in a failure to find the mesh file and write the output files unless the paths (in the json file) are explicitly provided as an input + # to the process. + # This is not an issue in the case of Snakemake, as the working directory doesn't automatically change between different rules. + input: + parameters = lambda wildcards: configuration_to_parameter_file[wildcards.configuration], + mdpa = f"{result_dir}/mesh_{{configuration}}.mdpa", + kratos_input_template = f"{tool}/input_template.json", + kratos_material_template = f"{tool}/StructuralMaterials_template.json", + script_create_kratos_input = f"{tool}/create_kratos_input.py", + script_run_kratos_simulation = f"{tool}/run_kratos_simulation.py", + output: + kratos_inputfile = f"{result_dir}/ProjectParameters_{{configuration}}.json", + kratos_materialfile = f"{result_dir}/MaterialParameters_{{configuration}}.json", + result_vtk = f"{result_dir}/{{configuration}}/Structure_0_1.vtk", + conda: + f"{tool}/environment_simulation.yml", + shell: + """ + python3 {input.script_create_kratos_input} \ + --input_parameter_file {input.parameters} \ + --input_mdpa_file {input.mdpa} \ + --input_kratos_input_template {input.kratos_input_template} \ + --input_material_template {input.kratos_material_template} \ + --output_kratos_inputfile {output.kratos_inputfile} \ + --output_kratos_materialfile {output.kratos_materialfile} + + python3 {input.script_run_kratos_simulation} \ + --input_parameter_file {input.parameters} \ + --input_kratos_inputfile {output.kratos_inputfile} \ + --input_kratos_materialfile {output.kratos_materialfile} + """ + +rule postprocess_kratos_results: + input: + parameters = lambda wildcards: configuration_to_parameter_file[wildcards.configuration], + result_vtk = f"{result_dir}/{{configuration}}/Structure_0_1.vtk", + script = f"{tool}/postprocess_results.py", + output: + zip = f"{result_dir}/solution_field_data_{{configuration}}.zip", + metrics = f"{result_dir}/solution_metrics_{{configuration}}.json", + conda: + f"{tool}/environment_simulation.yml", + shell: + """ + python3 {input.script} \ + --input_parameter_file {input.parameters} \ + --input_result_vtk {input.result_vtk} \ + --output_solution_file_zip {output.zip} \ + --output_metrics_file {output.metrics} + """ \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/StructuralMaterials_template.json b/examples/linear-elastic-plate-with-hole/kratos/kratos/StructuralMaterials_template.json new file mode 100644 index 0000000..4b51dd9 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/StructuralMaterials_template.json @@ -0,0 +1,18 @@ +{ + "properties": [ + { + "model_part_name": "Structure", + "properties_id": 1, + "Material": { + "constitutive_law": { + "name": "LinearElasticPlaneStress2DLaw" + }, + "Variables": { + "YOUNG_MODULUS": "{{YOUNG_MODULUS}}", + "POISSON_RATIO": "{{POISSON_RATIO}}" + }, + "Tables": {} + } + } + ] +} \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/create_kratos_input.py b/examples/linear-elastic-plate-with-hole/kratos/kratos/create_kratos_input.py new file mode 100644 index 0000000..bb6b0c8 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/create_kratos_input.py @@ -0,0 +1,141 @@ +import json +import os +from pint import UnitRegistry +from argparse import ArgumentParser +from pathlib import Path +import sys +# Ensure the parent directory is in the path to import AnalyticalSolution +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from analytical_solution import AnalyticalSolution + +def create_kratos_input( + parameter_file: str, + mdpa_file: str, + kratos_input_template_file: str, + kratos_material_template_file: str, + kratos_input_file: str, + kratos_material_file: str, +): + ureg = UnitRegistry() + with open(parameter_file) as f: + parameters = json.load(f) + + E = ( + ureg.Quantity( + parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"] + ) + .to_base_units() + .magnitude + ) + nu = ( + ureg.Quantity( + parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"] + ) + .to_base_units() + .magnitude + ) + radius = ( + ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + .to_base_units() + .magnitude + ) + L = ( + ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + .to_base_units() + .magnitude + ) + load = ( + ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"]) + .to_base_units() + .magnitude + ) + + analytical_solution = AnalyticalSolution( + E=E, + nu=nu, + radius=radius, + L=L, + load=load, + ) + + bc = analytical_solution.displacement_symbolic_str("X", "Y") + + with open(kratos_material_template_file) as f: + material_string = f.read() + + material_string = material_string.replace(r'"{{YOUNG_MODULUS}}"', str(E)) + material_string = material_string.replace(r'"{{POISSON_RATIO}}"', str(nu)) + + with open(kratos_material_file, "w") as f: + f.write(material_string) + + with open(kratos_input_template_file) as f: + project_parameters_string = f.read() + project_parameters_string = project_parameters_string.replace( + r"{{MESH_FILE}}", os.path.splitext(mdpa_file)[0] + ) + project_parameters_string = project_parameters_string.replace( + r"{{MATERIAL_FILE}}", kratos_material_file + ) + project_parameters_string = project_parameters_string.replace( + r"{{BOUNDARY_RIGHT_DISPLACEMENT_X}}", str(bc[0]) + ) + project_parameters_string = project_parameters_string.replace( + r"{{BOUNDARY_RIGHT_DISPLACEMENT_Y}}", str(bc[1]) + ) + project_parameters_string = project_parameters_string.replace( + r"{{BOUNDARY_TOP_DISPLACEMENT_X}}", str(bc[0]) + ) + project_parameters_string = project_parameters_string.replace( + r"{{BOUNDARY_TOP_DISPLACEMENT_Y}}", str(bc[1]) + ) + config = parameters["configuration"] + output_dir = os.path.join(os.path.dirname(os.path.abspath(kratos_input_file)), str(config)) + os.makedirs(output_dir, exist_ok=True) + project_parameters_string = project_parameters_string.replace(r"{{OUTPUT_PATH}}", output_dir) + + with open(kratos_input_file, "w") as f: + f.write(project_parameters_string) + +if __name__ == "__main__": + parser = ArgumentParser( + description="Create Kratos input and material files from templates and parameters." + ) + parser.add_argument( + "--input_parameter_file", + required=True, + help="JSON file containing simulation parameters (input)", + ) + parser.add_argument( + "--input_mdpa_file", required=True, help="Path to the MDPA mesh file (input)" + ) + parser.add_argument( + "--input_kratos_input_template", + required=True, + help="Path to the kratos input template file (input)", + ) + parser.add_argument( + "--input_material_template", + required=True, + help="Path to the kratos material template file (input)", + ) + parser.add_argument( + "--output_kratos_inputfile", + required=True, + help="Path to the kratos input file (output)", + ) + parser.add_argument( + "--output_kratos_materialfile", + required=True, + help="Path to the kratos material file (output)", + ) + args, _ = parser.parse_known_args() + + create_kratos_input( + parameter_file=args.input_parameter_file, + mdpa_file=args.input_mdpa_file, + kratos_input_template_file=args.input_kratos_input_template, + kratos_material_template_file=args.input_material_template, + kratos_input_file=args.output_kratos_inputfile, + kratos_material_file=args.output_kratos_materialfile, + ) diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/environment_simulation.yml b/examples/linear-elastic-plate-with-hole/kratos/kratos/environment_simulation.yml new file mode 100644 index 0000000..0bee2ef --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/environment_simulation.yml @@ -0,0 +1,15 @@ +name: kratos_simulation +# Environment file for kratos simulation scripts. Called by kratos tool workflow. + +channels: + - conda-forge +dependencies: + - python=3.10 + - meshio + - python-gmsh + - pint + - pyvista + - pip + - pip: + - KratosMultiphysics-all + - sympy diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/input_template.json b/examples/linear-elastic-plate-with-hole/kratos/kratos/input_template.json new file mode 100644 index 0000000..f8b680a --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/input_template.json @@ -0,0 +1,142 @@ +{ + "problem_data": { + "problem_name": "PlateWithHole", + "parallel_type": "OpenMP", + "start_time": 0.0, + "end_time": 1.0, + "echo_level": 0 + }, + "solver_settings": { + "solver_type": "Static", + "model_part_name": "Structure", + "echo_level": 1, + "domain_size": 2, + "analysis_type": "linear", + "model_import_settings": { + "input_type": "mdpa", + "input_filename": "{{MESH_FILE}}" + }, + "material_import_settings": { + "materials_filename": "{{MATERIAL_FILE}}" + }, + "time_stepping": { + "time_step": 1.0 + } + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "Parameters": { + "model_part_name": "Structure.boundary_left", + "variable_name": "DISPLACEMENT", + "constrained": [ + true, + false, + true + ], + "value": [ + 0.0, + 0.0, + 0.0 + ], + "interval": [ + 0.0, + "End" + ] + } + }, + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "Parameters": { + "model_part_name": "Structure.boundary_bottom", + "variable_name": "DISPLACEMENT", + "constrained": [ + false, + true, + true + ], + "value": [ + 0.0, + 0.0, + 0.0 + ], + "interval": [ + 0.0, + "End" + ] + } + }, + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "Parameters": { + "model_part_name": "Structure.boundary_right", + "variable_name": "DISPLACEMENT", + "constrained": [ + true, + true, + true + ], + "value": [ + "{{BOUNDARY_RIGHT_DISPLACEMENT_X}}", + "{{BOUNDARY_RIGHT_DISPLACEMENT_Y}}", + 0.0 + ], + "interval": [ + 0.0, + "End" + ] + } + }, + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "Parameters": { + "model_part_name": "Structure.boundary_top", + "variable_name": "DISPLACEMENT", + "constrained": [ + true, + true, + true + ], + "value": [ + "{{BOUNDARY_TOP_DISPLACEMENT_X}}", + "{{BOUNDARY_TOP_DISPLACEMENT_Y}}", + 0.0 + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "loads_process_list": [], + "list_other_processes": [] + }, + "output_processes": { + "vtk_output": [ + { + "python_module": "vtk_output_process", + "kratos_module": "KratosMultiphysics", + "Parameters": { + "model_part_name": "Structure", + "file_format": "binary", + "output_path": "{{OUTPUT_PATH}}", + "output_sub_model_parts": false, + "output_interval": 1, + "nodal_solution_step_data_variables": [ + "DISPLACEMENT" + ], + "gauss_point_variables_extrapolated_to_nodes": [ + "CAUCHY_STRESS_VECTOR", + "VON_MISES_STRESS" + ] + } + } + ] + } +} \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/msh_to_mdpa.py b/examples/linear-elastic-plate-with-hole/kratos/kratos/msh_to_mdpa.py new file mode 100644 index 0000000..78a135b --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/msh_to_mdpa.py @@ -0,0 +1,132 @@ +import json +import meshio +import re +import numpy as np +from pint import UnitRegistry +from argparse import ArgumentParser + +def msh_to_mdpa(parameter_file: str, mesh_file: str, mdpa_file: str): + """ + This function converts the GMSH mesh to a Kratos MDPA file format. + Due to limitations in the meshio conversion, several modifications are made to + the mdpa file: + - The element types are replaced with SmallDisplacementElement2D3N and SmallDisplacementElement2D6N + since meshio only converts to Triangle2D3 and Triangle2D6 which only describe the geometry but + not the finite elements. + - The Line2D elements are removed since they are not used in Kratos. + - The gmsh:dim_tags are removed since they are not used in Kratos. + - SubModelParts for the boundary conditions are created. + + At this point, we don't see a better way to do this conversion, so we use a lot of string manipulation. + """ + + ureg = UnitRegistry() + with open(parameter_file) as f: + parameters = json.load(f) + radius = ( + ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + .to_base_units() + .magnitude + ) + L = ( + ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + .to_base_units() + .magnitude + ) + + x0 = 0.0 + x1 = x0 + radius + x2 = x0 + L + y0 = 0.0 + y1 = y0 + radius + y2 = y0 + L + mesh = meshio.read(mesh_file) + + meshio.write(mdpa_file, mesh) + + with open(mdpa_file, "r") as f: + # replace all occurences of Triangle with SmallStrainElement + text = f.read() + + text = text.replace("Triangle2D3", "SmallDisplacementElement2D3N") + text = text.replace("Triangle2D6", "SmallDisplacementElement2D6N") + + text = re.sub(r"Begin\s+Elements\s+Line2D[\n\s\d]*End\s+Elements", "", text) + + mesh_tags = np.array( + re.findall( + r"Begin\s+NodalData\s+gmsh:dim_tags[\s\n]*(.*)End\s+NodalData\s+gmsh:dim_tags", + text, + flags=re.DOTALL, + )[0] + .replace("np.int64", "") + .replace("(", "") + .replace(")", "") + .split(), + dtype=np.int32, + ).reshape(-1, 3) + + text = re.sub( + r"Begin\s+NodalData\s+gmsh:dim_tags[\s\n]*(.*)End\s+NodalData\s+gmsh:dim_tags", + "", + text, + flags=re.DOTALL, + ) + + append = "\nBegin SubModelPart boundary_left\n" + append += " Begin SubModelPartNodes\n " + nodes = np.argwhere(np.isclose(mesh.points[:, 0], x0)).flatten() + 1 + append += "\n ".join(map(str, nodes)) + "\n" + append += " End SubModelPartNodes\n" + append += "End SubModelPart\n" + + text += append + + append = "\nBegin SubModelPart boundary_bottom\n" + append += " Begin SubModelPartNodes\n " + nodes = np.argwhere(np.isclose(mesh.points[:, 1], y0)).flatten() + 1 + append += "\n ".join(map(str, nodes)) + "\n" + append += " End SubModelPartNodes\n" + append += "End SubModelPart\n" + + text += append + + append = "\nBegin SubModelPart boundary_right\n" + append += " Begin SubModelPartNodes\n " + nodes = np.argwhere(np.isclose(mesh.points[:, 0], x2)).flatten() + 1 + append += "\n ".join(map(str, nodes)) + "\n" + append += " End SubModelPartNodes\n" + append += "End SubModelPart\n" + + text += append + + append = "\nBegin SubModelPart boundary_top\n" + append += " Begin SubModelPartNodes\n " + nodes = np.argwhere(np.isclose(mesh.points[:, 1], y2)).flatten() + 1 + append += "\n ".join(map(str, nodes)) + "\n" + append += " End SubModelPartNodes\n" + append += "End SubModelPart\n" + + text += append + with open(mdpa_file, "w") as f: + f.write(text) + +if __name__ == "__main__": + parser = ArgumentParser( + description="Convert GMSH mesh to Kratos MDPA format." + ) + parser.add_argument( + "--input_parameter_file", + required=True, + help="JSON file containing simulation parameters (input)", + ) + parser.add_argument( + "--input_mesh_file", required=True, help="Path to the mesh file (input)" + ) + parser.add_argument( + "--output_mdpa_file", + required=True, + help="Path to the MDPA file (output)", + ) + args, _ = parser.parse_known_args() + msh_to_mdpa(args.input_parameter_file, args.input_mesh_file, args.output_mdpa_file) diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/postprocess_results.py b/examples/linear-elastic-plate-with-hole/kratos/kratos/postprocess_results.py new file mode 100644 index 0000000..467fb42 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/postprocess_results.py @@ -0,0 +1,58 @@ +import json +import pyvista +from pathlib import Path +import zipfile +from argparse import ArgumentParser + +def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_file, output_solution_file_zip): + with open(input_parameter_file) as f: + parameters = json.load(f) + config = parameters["configuration"] + + mesh = pyvista.read(str(input_result_vtk)) + max_von_mises_stress = float(mesh["VON_MISES_STRESS"].max()) + print("Max Von Mises Stress:", max_von_mises_stress) + metrics = { + "max_von_mises_stress_nodes": max_von_mises_stress + } + with open(output_metrics_file, "w") as f: + json.dump(metrics, f, indent=4) + + files_to_store = [str(input_result_vtk)] + + with zipfile.ZipFile(output_solution_file_zip, "w") as zipf: + for filepath in files_to_store: + zipf.write(filepath, arcname=f"result_{config}.vtk") + +if __name__ == "__main__": + parser = ArgumentParser( + description="Postprocess Kratos results and write metrics and zipped solution." + ) + parser.add_argument( + "--input_parameter_file", + required=True, + help="JSON file containing simulation parameters (input)", + ) + parser.add_argument( + "--input_result_vtk", + required=True, + help="Path to the Kratos result VTK file (input)", + ) + parser.add_argument( + "--output_solution_file_zip", + required=True, + help="Path to the zipped solution files (output)", + ) + parser.add_argument( + "--output_metrics_file", + required=True, + help="Path to the output metrics JSON file (output)", + ) + args, _ = parser.parse_known_args() + + postprocess_results( + args.input_parameter_file, + args.input_result_vtk, + args.output_metrics_file, + args.output_solution_file_zip + ) diff --git a/examples/linear-elastic-plate-with-hole/kratos/kratos/run_kratos_simulation.py b/examples/linear-elastic-plate-with-hole/kratos/kratos/run_kratos_simulation.py new file mode 100644 index 0000000..5df8a0b --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/kratos/run_kratos_simulation.py @@ -0,0 +1,51 @@ +from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 +import json +import sys +from argparse import ArgumentParser + +import gmsh +import meshio +import re +from pint import UnitRegistry +import numpy as np +import os + +ureg = UnitRegistry() + +import KratosMultiphysics +from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis +import sys + +import pyvista +from pathlib import Path + + +if __name__ == "__main__": + parser = ArgumentParser( + description="Run FEniCS simulation for a plate with a hole.\n" + "Inputs: --input_parameter_file, --input_kratos_inputfile, --input_kratos_materialfile\n" + "Outputs: --output_solution_file_zip, --output_metrics_file" + ) + parser.add_argument( + "--input_parameter_file", + required=True, + help="JSON file containing simulation parameters (input)", + ) + parser.add_argument( + "--input_kratos_inputfile", + required=True, + help="Path to the kratos input file (input)", + ) + parser.add_argument( + "--input_kratos_materialfile", + required=True, + help="Path to the kratos material file (input)", + ) + args, _ = parser.parse_known_args() + + with open(args.input_kratos_inputfile, "r") as kratos_input: + parameters = KratosMultiphysics.Parameters(kratos_input.read()) + + model = KratosMultiphysics.Model() + simulation = StructuralMechanicsAnalysis(model, parameters) + simulation.Run() \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py b/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py new file mode 100644 index 0000000..8883164 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py @@ -0,0 +1,63 @@ +from pathlib import Path +import zipfile +import json +from benchmarking import benchmark + +""" +The script performs the following steps: + +1. Extracts the benchmark files from a zip archive (currently assuming that it is an RO-Crate of the benchmark). +2. Creates a benchmark object and registers the simulation tool (FEniCS) along with its version and URI. +3. Adds the simulation script and environment file to the benchmark object. +4. Iterates through the parameter configuration files, checks the "element-size" value, and if it meets the specified condition (>= 0.025), it generates and executes the Snakemake workflow for that configuration. + +The results of each run (and the files used by it) are stored in the directory with the configuration name. +""" + + +#################################################################################################### +#################################################################################################### +# Benchmark Extraction +#################################################################################################### +#################################################################################################### + +zipped_benchmark_dir = Path(__file__).resolve().parent.parent +unzipped_benchmark_dir = Path(__file__).resolve().parent + +with zipfile.ZipFile(zipped_benchmark_dir / "linear-elastic-plate-with-hole.zip", 'r') as zip_ref: + # Extract all files + zip_ref.extractall(unzipped_benchmark_dir) + + +#################################################################################################### +#################################################################################################### +# Creation of benchmark object and the addition of simulation tool scripts +#################################################################################################### +#################################################################################################### + +linear_elastic_problem = benchmark("linear-elastic-plate-with-hole", \ + benchmark_dir=unzipped_benchmark_dir, \ + benchmark_uri="https://portal.mardi4nfdi.de/wiki/Model:6775296") + +linear_elastic_problem.add_tool("kratos", \ + version="10.3.1", \ + uri="https://github.com/KratosMultiphysics/Kratos") + +""" linear_elastic_problem.add_tool_scripts( + simulation_script = unzipped_benchmark_dir / "run_fenics_simulation.py", + environment_file = unzipped_benchmark_dir / "environment_simulation.yml" +) + +#################################################################################################### +#################################################################################################### +# Conditional creation and execution of snakemake workflows based on parameter configurations +#################################################################################################### +#################################################################################################### + +for file in unzipped_benchmark_dir.glob("parameters_*.json"): + with open(file, "r") as f: + data = json.load(f) + if data.get("element-size").get("value") >= 0.025: + linear_elastic_problem.generate_workflow(file.name, data.get("configuration")) + linear_elastic_problem.run_workflow() + """ \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip new file mode 100644 index 0000000..c6d6172 Binary files /dev/null and b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip differ diff --git a/pyproject.toml b/pyproject.toml index a477fe6..3d52e4f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,16 @@ [build-system] -requires = ["hatchling"] -build-backend = "hatchling.build" +requires = ["setuptools>=42", "wheel"] +build-backend = "setuptools.build_meta" [project] -name = "meshhelper" -version = "0.1.0" \ No newline at end of file +name = "benchmarking" +version = "0.1.0" +description = "Repository for model validation and verification" +authors = [ + { name="Your Name", email="you@example.com" } +] +readme = "README.md" +requires-python = ">=3.8" + +[tool.setuptools.packages.find] +where = ["src"] \ No newline at end of file diff --git a/src/__init__.py b/src/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/benchmarking.py b/src/benchmarking.py new file mode 100644 index 0000000..e11fe5a --- /dev/null +++ b/src/benchmarking.py @@ -0,0 +1,205 @@ +from pathlib import Path +from typing import Optional +import shutil +import subprocess + + + + +class benchmark: + """ + Represents a benchmark instance and manages all related operations. + """ + def __init__(self, name: str, benchmark_dir: Path, benchmark_uri: Optional[str] = None): + """ + Initialize a Benchmark instance. + + Args: + name: The name of the benchmark + benchmark_dir: Optional path to the benchmark directory + benchmark_uri: Optional URI/URL for the benchmark + """ + self.name = name + self.benchmark_dir = benchmark_dir + self.benchmark_uri = benchmark_uri + self.mesh = {} + self.tool = {} + + + + + def add_tool(self, name: str, version: Optional[str] = None, uri: Optional[str] = None): + """ + Register a new tool with optional metadata. + + Args: + name: Name of the tool + version: Optional version string + uri: Optional URI/URL for tool documentation + """ + + #if not Path(name).exists(): + # raise FileNotFoundError(f"Tool directory not found: {name}") + + tool_info = { + "name": name, + "version": version, + "uri": uri if uri is not None else "NA", + } + self.tool.update(tool_info) + + + + + def add_tool_scripts(self, simulation_script: Path, environment_file: Path): + """ + Add a command or script for the tool. + + Args: + tool_name: Name of the tool to add script to + script_cmd: Command or script path + environment_file: Path to conda environment file + """ + + #if not (Path(self.tool["name"]) / simulation_script).exists(): + if not simulation_script.exists(): + raise FileNotFoundError(f"Simulation script not found: {simulation_script}") + + #if not (Path(self.tool["name"]) / environment_file).exists(): + if not environment_file.exists(): + raise FileNotFoundError(f"Environment file not found: {environment_file}") + + self.tool.update({"simulation_script": simulation_script, "environment_file": environment_file}) + + + + + def add_tool_workflow(self, workflow_file: Path): + """ + Add a Snakemake workflow file for the tool. + + Args: + workflow_file: Path to the Snakemake workflow file + """ + + if workflow_file.exists(): + raise FileNotFoundError(f"Workflow file not found: {workflow_file}") + + self.tool.update({"workflow_file": workflow_file}) + + + + + def generate_workflow(self, parameter_file: str, configuration: str): + """ + Generate the Snakemake workflow for the benchmark with the associated tool. + parameter_file: The name of the parameter file to be used in the workflow. + configuration: The name of the configuration to be used for naming the output directory and files. + """ + + ############################################################################### + #Read the simulation tool workflow template + ############################################################################### + + if "simulation_script" in self.tool and "workflow_file" not in self.tool: + # Load template from external file + tool_template_path = self.benchmark_dir / "tool_workflow_template.txt" + + with open(tool_template_path, 'r') as f: + tool_template = f.read() + + # Replace placeholders with actual values + tool_workflow_content = tool_template.replace("$SIMULATION_SCRIPT$", self.tool["simulation_script"].name) \ + .replace("$TOOL_ENVIRONMENT_FILE$", self.tool["environment_file"].name) + + ############################################################################### + #Read the simulation tool user-defined workflow + ############################################################################### + + elif "workflow_file" in self.tool and "simulation_script" not in self.tool: + #tool_workflow_path = Path(self.tool["name"]) / self.tool["workflow_file"] + with open(self.tool["workflow_file"], 'r') as f: + tool_workflow_content = f.read() + + else: + raise ValueError("Either tool scripts (the simulation script and the environment file) or tool workflow must be provided.") + + ############################################################################### + # Append the simulation tool workflow to the main workflow template + ############################################################################### + + main_template_path = self.benchmark_dir / "main_workflow_template.txt" + with open(main_template_path, 'r') as f: + main_template = f.read() + + main_workflow_content = main_template.replace("$TOOL_WORKFLOW$", tool_workflow_content) \ + .replace("$BENCHMARK_NAME$", self.name) \ + .replace("$BENCHMARK_URI$", "https://portal.mardi4nfdi.de/wiki/Model:6775296") + + + self.output_dir = self.benchmark_dir / "results" / configuration + self.output_dir.mkdir(parents=True, exist_ok=True) + + # Copy files from benchmark_dir to self.output_dir, excluding non-matching parameter files and workflow template files + for item in self.benchmark_dir.iterdir(): + if item.is_file(): + if item.name.startswith("parameters_"): + # Only copy the matching parameter file + if item.name == parameter_file: + item_copy = self.output_dir / "parameters.json" + shutil.copy(item, item_copy) + elif item.name not in [tool_template_path.name, main_template_path.name]: # Exclude template files + # Copy all non-parameter files + shutil.copy(item, self.output_dir / item.name) + + + with open(self.output_dir / "Snakefile", 'w') as f: + f.write(main_workflow_content) + + print(f"Snakefile generated successfully") + + + + + def run_workflow(self,): + """ + Run the generated Snakemake workflow. + """ + if not (self.output_dir / "Snakefile").exists(): + raise ValueError("Snakemake workflow file not found. Please generate the workflow first.") + + #Creates a directory to store the conda environments. The environments are shared across different parameter configurations. + #To avoid redundant creation of environments, this path will be passed to all snakemake files during execution. + + shared_env_dir = self.benchmark_dir / "conda_envs" + shared_env_dir.mkdir(parents=True, exist_ok=True) + + # Run the Snakemake workflow using subprocess + try: + subprocess.run(["snakemake", "--use-conda", "--force", "--cores", "all", "--conda-prefix", str(self.benchmark_dir / "conda_envs")], check=True, cwd=self.output_dir) + print("Workflow executed successfully.") + + #Running the reporter plugin. + #subprocess.run(["snakemake", "--use-conda", "--force", "--cores", "all", \ + # "--reporter", "metadata4ing", \ + # "--report-metadata4ing-paramscript", "../common/parameter_extractor.py",\ + # "--report-metadata4ing-config", "metadata4ing.config", \ + # "--report-metadata4ing-filename", "$SNAKEMAKE_RESULT_FILE"], check=True, cwd=self.output_dir) + + except subprocess.CalledProcessError as e: + print(f"Error occurred while running the workflow: {e}") + + + + + + + + + + + + + + +