diff --git a/.github/workflows/run-benchmark.yml b/.github/workflows/run-benchmark.yml index 546edb1..b4cec3a 100644 --- a/.github/workflows/run-benchmark.yml +++ b/.github/workflows/run-benchmark.yml @@ -54,12 +54,6 @@ jobs: --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 - - name: Archive Linear Elastic plate with a hole benchmark data for snakemake uses: actions/upload-artifact@v4 with: @@ -67,13 +61,6 @@ jobs: 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/ - process-artifacts: runs-on: ubuntu-latest needs: run-simulation diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..4b5a294 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,4 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:conda", + "python-envs.defaultPackageManager": "ms-python.python:conda" +} \ No newline at end of file diff --git a/benchmarks/common/nextflow.config b/benchmarks/common/nextflow.config deleted file mode 100644 index 36401ae..0000000 --- a/benchmarks/common/nextflow.config +++ /dev/null @@ -1,27 +0,0 @@ -conda { - enabled = true -} - -params.result_dir = "nextflow_results/${params.benchmark}" - -prov { - formats { - dag { - file = "${params.result_dir}/nf_prov_dag.html" - overwrite = true - } - legacy { - file = "${params.result_dir}/nf_prov_legacy.json" - overwrite = true - } - wrroc { - agent { - name = 'Firstname Lastname' - orcid = 'https://orcid.org/0000-0000-0000-0000' - } - file = "${params.result_dir}/ro-crate-metadata.json" - license = 'https://spdx.org/licenses/MIT' - overwrite = true - } - } -} 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/docs/benchmarks/linear elasticity/plate-with-hole.md b/docs/benchmarks/linear elasticity/plate-with-hole.md index fcc8861..cb37073 100644 --- a/docs/benchmarks/linear elasticity/plate-with-hole.md +++ b/docs/benchmarks/linear elasticity/plate-with-hole.md @@ -94,36 +94,40 @@ B(\boldsymbol u,\boldsymbol v) &= \int_{\Omega} \boldsymbol\varepsilon(\boldsymb $$ -In order to solve the weak formulation, the finite-element method (FEM) can be used. This method discretizes the domain $\Omega$ into so called finite elements that can for example be triangles or quadrilaterals in 2D. On these elements, ansatz functions are defined such that they are continous on the boundaries between elements. These functions form a basis for the solution space for an approximate solution $\boldsymbol{u}_h$ of the problem. +In order to solve the weak formulation, the finite-element method (FEM) can be used. This method discretizes the domain $\Omega$ into so called finite elements that can for example be triangles or quadrilaterals in 2D. On these elements, ansatz functions are defined such that they are continuous on the boundaries between elements. These functions form a basis for the solution space for an approximate solution $\boldsymbol{u}_h$ of the problem. -## Comparison of approximate solution with analytical solution +## Output metrics -The approxiamte solution and the analytical solution can be compared with the $L_2$ norm which is defined as +- The approximate solution is compared with the analytical solution using the $L_2$ norm which is defined as -$$ -\Vert \boldsymbol{u}\Vert_{L_2} = \sqrt{\int_\Omega \Vert\boldsymbol{u}(\boldsymbol{x})\Vert_2^2 \mathrm d \boldsymbol x} -$$ + $$ + \Vert \boldsymbol{f}\Vert_{L_2} = \sqrt{\int_\Omega \left|\boldsymbol{f}(\boldsymbol{x})\right|^2 \mathrm d \boldsymbol x} + $$ -and the error in the $L_2$ norm + This norm is computed for the error between the finite element solution and the analytical solution of displacements i.e., $\Vert \boldsymbol{u}-\boldsymbol{u}_h\Vert_{L_2}$. -$$ -e_{L_2} = \Vert \boldsymbol{u}-\boldsymbol{u}_h\Vert_{L_2}. -$$ +- The maximum displacement error is computed at the nodes of the finite element mesh with respect to the analytical solution. -Alternatively, the sup norm can be used which is defined as + $$ + e_{\max} + = + \max_{i \in \mathcal{N}} + \left\| + \mathbf{u}(\mathbf{x}_i) + - + \mathbf{u}_h(\mathbf{x}_i) + \right\| + $$ -$$ -\Vert \boldsymbol{u}\Vert_{\inf} = \sup_{\boldsymbol x} \Vert \boldsymbol{u}(\boldsymbol{x}) \Vert -$$ + where $\mathcal{N}$ denotes the set of nodes of the finite element mesh. -and the error in the sup norm +- Max Von-Mises stress -$$ -e_{\inf} = \Vert \boldsymbol{u}-\boldsymbol{u}_h\Vert_{\inf}. -$$ +- Displacement at the top-right corner of the plate. +- The reaction force at the left boundary. -With these metrices, we can perform a convergence analysis for different approximations $\boldsymbol{u}_h$ which differ in the element size $h$. Plotting the error over the used element-size in a log-log plot lets us determine the convergence order of the approximation. +- The number of Degrees of Freedom in the finite element mesh. ## Table of parameters diff --git a/examples/linear-elastic-plate-with-hole/fenics/Snakefile b/examples/linear-elastic-plate-with-hole/fenics/Snakefile new file mode 100644 index 0000000..b463a9f --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/fenics/Snakefile @@ -0,0 +1,33 @@ +import json + +rule all: + input: + "solution_metrics.json", + "solution_field_data.zip" + +rule create_mesh: + input: + script = "create_mesh.py", + parameters = "parameters.json", + output: + mesh = "mesh.msh", + conda: "environment_mesh.yml" + shell: + """ + python3 {input.script} --input_parameter_file {input.parameters} --output_mesh_file {output.mesh} + """ + +rule run_simulation: + input: + script = "run_simulation.py", + parameters = "parameters.json", + mesh = "mesh.msh", + output: + zip = "solution_field_data.zip", + metrics = "solution_metrics.json", + conda: + "environment_simulation.yml" + shell: + """ + python3 {input.script} --input_parameter_file {input.parameters} --input_mesh_file {input.mesh} --output_solution_file_zip {output.zip} --output_metrics_file {output.metrics} + """ diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py index 22d1f5f..852fec5 100644 --- a/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py +++ b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py @@ -34,6 +34,16 @@ shared_env_dir = root_unzipped_benchmark_dir / "conda_envs" shared_env_dir.mkdir(parents=True, exist_ok=True) +#################################################################################################### +#################################################################################################### +# Simulation tool metadata (to be included in the RO-Crate) +#################################################################################################### +#################################################################################################### + +tool_name = "fenics" +tool_uri = "https://github.com/FEniCS/dolfinx" +tool_version = "0.9.0" + #################################################################################################### #################################################################################################### # Conditional execution of parameter configurations @@ -43,7 +53,7 @@ for file in root_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: + if data.get("element_size[m]") >= 0.025: # Create output directory for the configuration output_dir = root_unzipped_benchmark_dir / "results" / data.get("configuration") diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_simulation.py b/examples/linear-elastic-plate-with-hole/fenics/run_simulation.py index 1aec6e5..eee0b1d 100644 --- a/examples/linear-elastic-plate-with-hole/fenics/run_simulation.py +++ b/examples/linear-elastic-plate-with-hole/fenics/run_simulation.py @@ -10,6 +10,7 @@ from dolfinx.fem.petsc import LinearProblem from mpi4py import MPI from pint import UnitRegistry +from ufl.algorithms import estimate_total_polynomial_degree # Add parent directory to sys.path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) @@ -29,7 +30,7 @@ def run_fenics_simulation( gdim=2, ) - V = df.fem.functionspace(mesh, ("CG", parameters["element-degree"], (2,))) + V = df.fem.functionspace(mesh, ("CG", parameters["isoparametric_element_degree"], (2,))) tags_left = facet_tags.find(1) tags_bottom = facet_tags.find(2) @@ -45,47 +46,30 @@ def run_fenics_simulation( E = ( ureg.Quantity( - parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"] + parameters["youngs_modulus[Pa]"], "Pa" ) .to_base_units() .magnitude ) nu = ( ureg.Quantity( - parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"] + parameters["poissons_ratio"], "" ) .to_base_units() .magnitude ) radius = ( - ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + ureg.Quantity(parameters["radius[m]"], "m") .to_base_units() .magnitude ) L = ( - ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + ureg.Quantity(parameters["length[m]"], "m") .to_base_units() .magnitude ) load = ( - ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"]) - .to_base_units() - .magnitude - ) - displacement_evaluation_point = parameters["displacement-evaluation-point"] - displacement_evaluation_x = ( - ureg.Quantity( - displacement_evaluation_point["x"]["value"], - displacement_evaluation_point["x"]["unit"], - ) - .to_base_units() - .magnitude - ) - displacement_evaluation_y = ( - ureg.Quantity( - displacement_evaluation_point["y"]["value"], - displacement_evaluation_point["y"]["unit"], - ) + ureg.Quantity(parameters["load[Pa]"], "Pa") .to_base_units() .magnitude ) @@ -115,10 +99,6 @@ def as_tensor(v): dx = ufl.Measure( "dx", - metadata={ - "quadrature_degree": parameters["quadrature-degree"], - "quadrature_scheme": parameters["quadrature-rule"], - }, ) ds = ufl.Measure( "ds", @@ -173,6 +153,7 @@ def traction_top_expr(x: np.ndarray) -> np.ndarray: reaction_left_y_local = df.fem.assemble_scalar(df.fem.form(traction[1] * ds(1))) reaction_left_x = MPI.COMM_WORLD.allreduce(reaction_left_x_local, op=MPI.SUM) reaction_left_y = MPI.COMM_WORLD.allreduce(reaction_left_y_local, op=MPI.SUM) + num_dofs = V.dofmap.index_map.size_global * V.dofmap.index_map_bs # Compute L2 error norm between FE displacement and analytical displacement. @@ -190,9 +171,17 @@ def analytical_displacement_expr(x: np.ndarray) -> np.ndarray: l2_error_sq_global = MPI.COMM_WORLD.allreduce(l2_error_sq_local, op=MPI.SUM) l2_error_displacement = np.sqrt(l2_error_sq_global) + # Compute max nodal displacement error magnitude (global across MPI) + block_size = V.dofmap.index_map_bs + nodal_error = (u.x.array - u_analytical.x.array).reshape(-1, block_size) + max_displacement_error_nodes_local = np.max(np.linalg.norm(nodal_error, axis=1)) + max_displacement_error_nodes = MPI.COMM_WORLD.allreduce( + max_displacement_error_nodes_local, op=MPI.MAX + ) + # Evaluate displacement at the specified evaluation point displacement_eval_point = np.array( - [[displacement_evaluation_x, displacement_evaluation_y, 0.0]], + [[1.0, 1.0, 0.0]], dtype=np.float64, ) tree = df.geometry.bb_tree(mesh, mesh.topology.dim) @@ -202,12 +191,13 @@ def analytical_displacement_expr(x: np.ndarray) -> np.ndarray: colliding_cells = df.geometry.compute_colliding_cells( mesh, cell_candidates, displacement_eval_point ) - local_displacement_x = None + local_displacement = None if len(colliding_cells.links(0)) > 0: cell = colliding_cells.links(0)[0] - local_displacement_x = u.eval( + # u.eval returns a 2D array: shape (num_points, value_size) + local_displacement = u.eval( displacement_eval_point, np.array([cell], dtype=np.int32) - )[0] + ).tolist() # [ux, uy] def project( @@ -238,10 +228,10 @@ def project( return uh plot_space_stress = df.fem.functionspace( - mesh, ("DG", parameters["element-degree"] - 1, (2, 2)) + mesh, ("DG", parameters["isoparametric_element_degree"] - 1, (2, 2)) ) plot_space_mises = df.fem.functionspace( - mesh, ("DG", parameters["element-degree"] - 1, (1,)) + mesh, ("DG", parameters["isoparametric_element_degree"] - 1, (1,)) ) stress_nodes_red = project(sigma(u), plot_space_stress, dx) stress_nodes_red.name = "stress" @@ -284,14 +274,13 @@ def mises_stress(u): ) 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"], + scheme="default", + degree=estimate_total_polynomial_degree(u_), ) Q_mises = df.fem.functionspace(mesh, quad_element) @@ -302,31 +291,32 @@ def mises_stress(u): np.max(mises_qp.x.array), op=MPI.MAX ) - displacement_x_at_evaluation_point = None + displacement_at_evaluation_point = None if MPI.COMM_WORLD.rank == 0: - displacement_x_candidates = ( - MPI.COMM_WORLD.gather(local_displacement_x, root=0) or [] + displacement_candidates = ( + MPI.COMM_WORLD.gather(local_displacement, root=0) or [] ) - for value in displacement_x_candidates: + for value in displacement_candidates: if value is not None: - displacement_x_at_evaluation_point = value + displacement_at_evaluation_point = value break - if displacement_x_at_evaluation_point is None: + if displacement_at_evaluation_point is None: raise ValueError( "Could not evaluate displacement at the configured evaluation point." ) else: - MPI.COMM_WORLD.gather(local_displacement_x, root=0) + MPI.COMM_WORLD.gather(local_displacement, root=0) # Save metrics metrics = { - "max_von_mises_stress_nodes": max_mises_stress_nodes, - "max_von_mises_stress_gauss_points": max_mises_stress_gauss_points, - "l2_error_displacement": l2_error_displacement, - "reaction_force_left_boundary_x": reaction_left_x, - "reaction_force_left_boundary_y": reaction_left_y, - f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})": displacement_x_at_evaluation_point, + "number_of_dofs[-]": num_dofs, + "max_von_mises_stress[Pa]": max_mises_stress_gauss_points, + "l2_error_displacement[m]": l2_error_displacement, + "max_displacement_error[m]": max_displacement_error_nodes, + "reaction_force_left_boundary_x[N]": reaction_left_x, + "reaction_force_left_boundary_y[N]": reaction_left_y, + "displacement_top_right_corner[m]": displacement_at_evaluation_point, # [ux, uy] } if MPI.COMM_WORLD.rank == 0: diff --git a/examples/linear-elastic-plate-with-hole/kratos/Snakefile_kratos.smk b/examples/linear-elastic-plate-with-hole/kratos/Snakefile.smk similarity index 88% rename from examples/linear-elastic-plate-with-hole/kratos/Snakefile_kratos.smk rename to examples/linear-elastic-plate-with-hole/kratos/Snakefile.smk index 79f1170..0940cc6 100644 --- a/examples/linear-elastic-plate-with-hole/kratos/Snakefile_kratos.smk +++ b/examples/linear-elastic-plate-with-hole/kratos/Snakefile.smk @@ -9,6 +9,18 @@ rule all: "solution_metrics.json", "solution_field_data.zip" +rule create_mesh: + input: + script = "create_mesh.py", + parameters = "parameters.json", + output: + mesh = "mesh.msh", + conda: "environment_mesh.yml" + shell: + """ + python3 {input.script} --input_parameter_file {input.parameters} --output_mesh_file {output.mesh} + """ + rule mesh_to_mdpa: input: parameters = "parameters.json", diff --git a/examples/linear-elastic-plate-with-hole/kratos/create_kratos_input.py b/examples/linear-elastic-plate-with-hole/kratos/create_kratos_input.py index f13c177..f9d727d 100644 --- a/examples/linear-elastic-plate-with-hole/kratos/create_kratos_input.py +++ b/examples/linear-elastic-plate-with-hole/kratos/create_kratos_input.py @@ -37,30 +37,30 @@ def create_kratos_input( E = ( ureg.Quantity( - parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"] + parameters["youngs_modulus[Pa]"], "Pa" ) .to_base_units() .magnitude ) nu = ( ureg.Quantity( - parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"] + parameters["poissons_ratio"], "" ) .to_base_units() .magnitude ) radius = ( - ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + ureg.Quantity(parameters["radius[m]"], "m") .to_base_units() .magnitude ) L = ( - ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + ureg.Quantity(parameters["length[m]"], "m") .to_base_units() .magnitude ) load = ( - ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"]) + ureg.Quantity(parameters["load[Pa]"], "Pa") .to_base_units() .magnitude ) diff --git a/examples/linear-elastic-plate-with-hole/kratos/msh_to_mdpa.py b/examples/linear-elastic-plate-with-hole/kratos/msh_to_mdpa.py index b07b965..f82c1a3 100644 --- a/examples/linear-elastic-plate-with-hole/kratos/msh_to_mdpa.py +++ b/examples/linear-elastic-plate-with-hole/kratos/msh_to_mdpa.py @@ -24,12 +24,12 @@ def msh_to_mdpa(parameter_file: str, mesh_file: str, mdpa_file: str): with open(parameter_file) as f: parameters = json.load(f) radius = ( - ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + ureg.Quantity(parameters["radius[m]"], "m") .to_base_units() .magnitude ) L = ( - ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + ureg.Quantity(parameters["length[m]"], "m") .to_base_units() .magnitude ) diff --git a/examples/linear-elastic-plate-with-hole/kratos/postprocess_results.py b/examples/linear-elastic-plate-with-hole/kratos/postprocess_results.py index 552e310..914aa00 100644 --- a/examples/linear-elastic-plate-with-hole/kratos/postprocess_results.py +++ b/examples/linear-elastic-plate-with-hole/kratos/postprocess_results.py @@ -24,48 +24,30 @@ def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_f E = ( ureg.Quantity( - parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"] + parameters["youngs_modulus[Pa]"], "Pa" ) .to_base_units() .magnitude ) nu = ( ureg.Quantity( - parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"] + parameters["poissons_ratio"], "" ) .to_base_units() .magnitude ) radius = ( - ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + ureg.Quantity(parameters["radius[m]"], "m") .to_base_units() .magnitude ) L = ( - ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + ureg.Quantity(parameters["length[m]"], "m") .to_base_units() .magnitude ) load = ( - ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"]) - .to_base_units() - .magnitude - ) - - displacement_evaluation_point = parameters["displacement-evaluation-point"] - displacement_evaluation_x = ( - ureg.Quantity( - displacement_evaluation_point["x"]["value"], - displacement_evaluation_point["x"]["unit"], - ) - .to_base_units() - .magnitude - ) - displacement_evaluation_y = ( - ureg.Quantity( - displacement_evaluation_point["y"]["value"], - displacement_evaluation_point["y"]["unit"], - ) + ureg.Quantity(parameters["load[Pa]"], "Pa") .to_base_units() .magnitude ) @@ -77,10 +59,8 @@ def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_f L=L, load=load, ) - # Compute maximum von Mises stress at nodes and Gauss points. - max_von_mises_stress_nodes = float(np.max(mesh.point_data["VON_MISES_STRESS"])) - - max_von_mises_stress_gauss_points = max_von_mises_stress_nodes + # Compute maximum von Mises stress at Gauss points. + max_von_mises_stress_gauss_points = 0 for key, values in mesh.cell_data.items(): if "VON_MISES_STRESS" in key: max_von_mises_stress_gauss_points = float(np.max(values)) @@ -112,24 +92,33 @@ def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_f reaction_force_left_boundary_x = float(np.sum(reaction[left_boundary_mask, 0])) reaction_force_left_boundary_y = float(np.sum(reaction[left_boundary_mask, 1])) + # Compute displacement at the top-right corner probe_points = pyvista.PolyData( - np.array([[displacement_evaluation_x, displacement_evaluation_y, 0.0]], dtype=float) + np.array([[1.0, 1.0, 0.0]], dtype=float) ) sampled = probe_points.sample(mesh) displacement_sampled = sampled.point_data.get("DISPLACEMENT") if displacement_sampled is None: - closest_id = mesh.find_closest_point([displacement_evaluation_x, displacement_evaluation_y, 0.0]) - displacement_x_at_evaluation_point = float(displacement[closest_id, 0]) + closest_id = mesh.find_closest_point([1.0, 1.0, 0.0]) + displacement_at_evaluation_point = [float(displacement[closest_id, 0]), float(displacement[closest_id, 1])] else: - displacement_x_at_evaluation_point = float(displacement_sampled[0, 0]) + displacement_at_evaluation_point = [float(displacement_sampled[0, 0]), float(displacement_sampled[0, 1])] + + # Compute nodal displacement error (Euclidean norm of error vector at each node) + nodal_displacement_error = np.linalg.norm(displacement - u_ref, axis=1) + max_displacement_error_nodes = float(np.max(nodal_displacement_error)) + + # Compute the number of dofs + num_dofs = int(mesh.n_points * 2) metrics = { - "max_von_mises_stress_nodes": max_von_mises_stress_nodes, - "max_von_mises_stress_gauss_points": max_von_mises_stress_gauss_points, - "l2_error_displacement": l2_error_displacement, - "reaction_force_left_boundary_x": reaction_force_left_boundary_x, - "reaction_force_left_boundary_y": reaction_force_left_boundary_y, - f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})": displacement_x_at_evaluation_point, + "number_of_dofs[-]": num_dofs, + "max_von_mises_stress[Pa]": max_von_mises_stress_gauss_points, + "l2_error_displacement[m]": l2_error_displacement, + "max_displacement_error[m]": max_displacement_error_nodes, + "reaction_force_left_boundary_x[N]": reaction_force_left_boundary_x, + "reaction_force_left_boundary_y[N]": reaction_force_left_boundary_y, + "displacement_top_right_corner[m]": displacement_at_evaluation_point, # [ux, uy] } with open(output_metrics_file, "w") as f: json.dump(metrics, f, indent=4) diff --git a/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py b/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py index 968d443..f417bb1 100644 --- a/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py +++ b/examples/linear-elastic-plate-with-hole/kratos/run_benchmark.py @@ -34,6 +34,16 @@ shared_env_dir = root_unzipped_benchmark_dir / "conda_envs" shared_env_dir.mkdir(parents=True, exist_ok=True) +#################################################################################################### +#################################################################################################### +# Simulation tool metadata (to be included in the RO-Crate) +#################################################################################################### +#################################################################################################### + +tool_name = "kratos" +tool_uri = "https://github.com/KratosMultiphysics/Kratos" +tool_version = "10.3.1" + #################################################################################################### #################################################################################################### # Conditional execution of parameter configurations @@ -43,7 +53,7 @@ for file in root_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: + if data.get("element_size[m]") >= 0.025: # Create output directory for the configuration output_dir = root_unzipped_benchmark_dir / "results" / data.get("configuration") @@ -60,11 +70,8 @@ continue else: shutil.copy(item, output_dir / item.name) - - # Run the Snakemake workflow from the benchmark to create the mesh for the configuration - subprocess.run(["snakemake", "--snakefile", "Snakefile", "--use-conda", "--force", "--cores", "all", "--conda-prefix", str(shared_env_dir), "--until", "create_mesh"], check=True, cwd=output_dir) - - # Run the Snakemake workflow to run the simulation and postprocess results for the configuration - subprocess.run(["snakemake", "--snakefile", "Snakefile_kratos.smk", "--use-conda", "--force", "--cores", "all", "--conda-prefix", str(shared_env_dir)], check=True, cwd=output_dir) + + # Run the Snakemake workflow for the configuration + subprocess.run(["snakemake", "--snakefile", "Snakefile.smk", "--use-conda", "--force", "--cores", "all", "--conda-prefix", str(shared_env_dir)], check=True, cwd=output_dir) print("Workflow executed successfully.") \ 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 index 0721633..080c1da 100644 Binary files a/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip and b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip differ