From 12643ed3c88b92031a28199721c2c24713930c1d Mon Sep 17 00:00:00 2001 From: IAntonau Date: Fri, 16 Aug 2024 11:59:40 +0200 Subject: [PATCH] [SiApp] p-norm damage resp (#12373) * save * adding p-norm to sensor sum * test case * revert files * apply changes suggested by Suneth * Update applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp Co-authored-by: Suneth Warnakulasuriya * save * adding testing * reverting ref values --------- Co-authored-by: IAntonau Co-authored-by: Suneth Warnakulasuriya --- .../add_custom_sensors_to_python.cpp | 2 +- ...measurement_residual_response_function.cpp | 30 +++- .../measurement_residual_response_function.h | 4 +- .../controls/material_properties_control.py | 2 +- .../responses/damage_detection_response.py | 2 +- .../system_identification_static_analysis.py | 5 +- .../adjoint_material_parameters.json | 3 + .../adjoint_project_parameters.json | 77 ++++++++++ .../optimization_parameters.json | 133 ++++++++++++++++++ .../sensor_data.json | 109 ++++++++++++++ ...stem_identification_p_norm_summary_ref.csv | 20 +++ .../tests/test_system_identification.py | 19 +++ 12 files changed, 396 insertions(+), 10 deletions(-) create mode 100644 applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_material_parameters.json create mode 100644 applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_project_parameters.json create mode 100644 applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/optimization_parameters.json create mode 100644 applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/sensor_data.json create mode 100644 applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm_summary_ref.csv diff --git a/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp b/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp index 3c4ffa27b90..866e9ad7110 100644 --- a/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp +++ b/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp @@ -99,7 +99,7 @@ void AddCustomSensorsToPython(pybind11::module& m) ; py::class_(sensor_module, "MeasurementResidualResponseFunction") - .def(py::init<>()) + .def(py::init(), py::arg("p_coefficient")) .def("AddSensor", &MeasurementResidualResponseFunction::AddSensor, py::arg("sensor")) .def("Clear", &MeasurementResidualResponseFunction::Clear) .def("GetSensorsList", &MeasurementResidualResponseFunction::GetSensorsList) diff --git a/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.cpp b/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.cpp index eec5fc2645c..fbdac3f7ee2 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.cpp +++ b/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.cpp @@ -6,7 +6,8 @@ // // License: SystemIdentificationApplication/license.txt // -// Main authors: Suneth Warnakulasuriya +// Main authors: Suneth Warnakulasuriya, +// Ihar Antonau // // System includes @@ -87,7 +88,8 @@ struct PartialSensitivity } // namespace MeasurementResidualResponseFunctionUtilities -MeasurementResidualResponseFunction::MeasurementResidualResponseFunction() +MeasurementResidualResponseFunction::MeasurementResidualResponseFunction(const double PCoefficient) + : mPCoefficient(PCoefficient) { mResponseGradientList.resize(ParallelUtilities::GetNumThreads()); } @@ -147,10 +149,13 @@ double MeasurementResidualResponseFunction::CalculateValue(ModelPart& rModelPart double value = 0.0; for (auto& p_sensor : mpSensorsList) { const double sensor_value = p_sensor->CalculateValue(rModelPart); - value += p_sensor->GetWeight() * std::pow(sensor_value - p_sensor->GetValue(SENSOR_MEASURED_VALUE), 2) * 0.5; p_sensor->SetSensorValue(sensor_value); + const double current_sensor_error_square = std::pow(sensor_value - p_sensor->GetValue(SENSOR_MEASURED_VALUE), 2) * 0.5; + p_sensor->SetValue(SENSOR_ERROR, current_sensor_error_square); + value += std::pow(p_sensor->GetWeight() * current_sensor_error_square, mPCoefficient); } - return value; + return std::pow(value, 1 / mPCoefficient); + KRATOS_CATCH(""); } @@ -170,11 +175,26 @@ void MeasurementResidualResponseFunction::CalculateDerivative( rResponseGradient.clear(); auto& local_sensor_response_gradient = mResponseGradientList[OpenMPUtils::ThisThread()]; + double temp = 0.0; + for (auto& p_sensor : mpSensorsList) { + temp += ( std::pow( p_sensor->GetValue(SENSOR_ERROR) * 0.5 * p_sensor->GetWeight(), mPCoefficient ) ); + } + const double c1 = 1 / mPCoefficient * std::pow( temp, 1/mPCoefficient - 1 ); + + temp = 0.0; + for (auto& p_sensor : mpSensorsList) { + temp += std::pow( p_sensor->GetWeight() * 0.5 * p_sensor->GetValue(SENSOR_ERROR), mPCoefficient - 1 ); + } + + const double c2 = mPCoefficient * temp; + for (auto& p_sensor : mpSensorsList) { TCalculationType::Calculate(*p_sensor, local_sensor_response_gradient, rResidualGradient, rArgs...); - noalias(rResponseGradient) += local_sensor_response_gradient * (p_sensor->GetWeight() * (p_sensor->GetSensorValue() - p_sensor->GetValue(SENSOR_MEASURED_VALUE))); + const double error = std::sqrt( p_sensor->GetValue(SENSOR_ERROR) ); + noalias(rResponseGradient) += c1 * c2 * error * local_sensor_response_gradient; } + KRATOS_CATCH(""); } diff --git a/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.h b/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.h index 2569bff1189..980964199c7 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.h +++ b/applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.h @@ -50,7 +50,7 @@ class KRATOS_API(DIGITAL_TWIN_APPLICATION) MeasurementResidualResponseFunction : ///@{ /// Constructor. - MeasurementResidualResponseFunction(); + MeasurementResidualResponseFunction(const double PCoeficient); /// Destructor. ~MeasurementResidualResponseFunction() override = default; @@ -155,6 +155,8 @@ class KRATOS_API(DIGITAL_TWIN_APPLICATION) MeasurementResidualResponseFunction : ///@name Member Variables ///@{ + double mPCoefficient; + std::vector mpSensorsList; std::vector mResponseGradientList; diff --git a/applications/SystemIdentificationApplication/python_scripts/controls/material_properties_control.py b/applications/SystemIdentificationApplication/python_scripts/controls/material_properties_control.py index f521c8e95c2..52163cc2da6 100644 --- a/applications/SystemIdentificationApplication/python_scripts/controls/material_properties_control.py +++ b/applications/SystemIdentificationApplication/python_scripts/controls/material_properties_control.py @@ -199,4 +199,4 @@ def _UpdateAndOutputFields(self, update: ContainerExpressionTypes) -> None: un_buffered_data.SetValue(f"{self.controlled_physical_variable.Name()}_physical_phi_derivative", self.physical_phi_derivative_field.Clone(), overwrite=True) def __str__(self) -> str: - return f"Control [type = {self.__class__.__name__}, name = {self.GetName()}, model part name = {self.adjoint_model_part_operation.GetModelPartFullName()}, control variable = {self.controlled_physical_variable.Name()}" + return f"Control [type = {self.__class__.__name__}, name = {self.GetName()}, model part name = {self.adjoint_model_part_operation.GetModelPartFullName()}, control variable = {self.controlled_physical_variable.Name()}" \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/python_scripts/responses/damage_detection_response.py b/applications/SystemIdentificationApplication/python_scripts/responses/damage_detection_response.py index e7f282ca529..f88de059aa8 100644 --- a/applications/SystemIdentificationApplication/python_scripts/responses/damage_detection_response.py +++ b/applications/SystemIdentificationApplication/python_scripts/responses/damage_detection_response.py @@ -92,7 +92,7 @@ def Initialize(self) -> None: def Check(self) -> None: pass - + def Finalize(self) -> None: self.adjoint_analysis.Finalize() diff --git a/applications/SystemIdentificationApplication/python_scripts/sensor_sensitivity_solvers/system_identification_static_analysis.py b/applications/SystemIdentificationApplication/python_scripts/sensor_sensitivity_solvers/system_identification_static_analysis.py index c8911ae3a14..1fccbfd03e2 100644 --- a/applications/SystemIdentificationApplication/python_scripts/sensor_sensitivity_solvers/system_identification_static_analysis.py +++ b/applications/SystemIdentificationApplication/python_scripts/sensor_sensitivity_solvers/system_identification_static_analysis.py @@ -22,6 +22,7 @@ def Initialize(self): "perturbation_size" : 1e-8, "adapt_perturbation_size" : true, "list_of_sensors" : [], + "p_coefficient" : 1, "output_settings" : { "output_sensor_sensitivity_fields": false, "output_folder" : "Optimization_Results/sensor_sensitivity_fields" @@ -36,7 +37,9 @@ def Initialize(self): model_part.ProcessInfo[KratosSI.ADAPT_PERTURBATION_SIZE] = sensor_settings["adapt_perturbation_size"].GetBool() self.listof_sensors = GetSensors(model_part, sensor_settings["list_of_sensors"].values()) - self.measurement_residual_response_function = KratosSI.Sensors.MeasurementResidualResponseFunction() + p_coefficient = sensor_settings["p_coefficient"].GetDouble() + self.measurement_residual_response_function = KratosSI.Sensors.MeasurementResidualResponseFunction(p_coefficient) + for sensor in self.listof_sensors: sensor.SetValue(KratosSI.SENSOR_MEASURED_VALUE, 0.0) self.measurement_residual_response_function.AddSensor(sensor) diff --git a/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_material_parameters.json b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_material_parameters.json new file mode 100644 index 00000000000..85b9c6d3473 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_material_parameters.json @@ -0,0 +1,3 @@ +{ + "properties": [] +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_project_parameters.json b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_project_parameters.json new file mode 100644 index 00000000000..ac2fd0dd756 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/adjoint_project_parameters.json @@ -0,0 +1,77 @@ +{ + "problem_data": { + "problem_name": "benchmark_2", + "parallel_type": "OpenMP", + "start_time": 0, + "end_time": 1, + "echo_level": 0 + }, + "sensor_settings": { + "perturbation_size": 1e-8, + "adapt_perturbation_size": true, + "p_coefficient": 4, + "@include_json": "auxiliary_files/system_identification/sensor_data.json" + }, + "solver_settings": { + "solver_type": "adjoint_static", + "analysis_type": "linear", + "model_part_name": "AdjointStructure", + "domain_size": 3, + "time_stepping": { + "time_step": 1.0 + }, + "compute_reactions": false, + "move_mesh_flag": false, + "sensitivity_settings": { + "sensitivity_model_part_name": "all_nodes_elements_model_part", + "element_data_value_sensitivity_variables": [ + "YOUNG_MODULUS" + ], + "build_mode": "static" + }, + "echo_level": 0, + "rotation_dofs": true, + "model_import_settings": { + "input_type": "use_input_model_part" + }, + "material_import_settings": { + "materials_filename": "auxiliary_files/system_identification/adjoint_material_parameters.json" + }, + "linear_solver_settings": { + "solver_type": "amgcl" + } + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "help": "This process fixes the selected components of a given vector variable", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "mesh_id": 0, + "model_part_name": "AdjointStructure.BoundPtsType7", + "variable_name": "ADJOINT_DISPLACEMENT", + "value": [ + 0.0, + 0.0, + 0.0 + ], + "constrained": [ + true, + true, + true + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "loads_process_list": [], + "list_other_processes": [] + }, + "output_processes": {}, + "analysis_stage": "KratosMultiphysics.SystemIdentificationApplication.sensor_sensitivity_solvers.sensor_sensitivity_static_analysis" +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/optimization_parameters.json b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/optimization_parameters.json new file mode 100644 index 00000000000..b2903ff8c6a --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/optimization_parameters.json @@ -0,0 +1,133 @@ +{ + "problem_data": { + "parallel_type": "OpenMP", + "echo_level": 0 + }, + "model_parts": [ + { + "type": "mdpa_model_part_controller", + "settings": { + "model_part_name": "AdjointStructure", + "domain_size": 3, + "input_filename": "auxiliary_files/structure" + } + }, + { + "type": "connectivity_preserving_model_part_controller", + "settings": { + "transformation_settings": [ + { + "source_model_part_name": "AdjointStructure", + "destination_model_part_name": "Structure", + "destination_element_name": "ShellThinElement3D3N", + "destination_condition_name": "LineLoadCondition3D2N" + } + ] + } + } + ], + "analyses": [ + { + "name": "Structure_static", + "type": "kratos_analysis_execution_policy", + "settings": { + "model_part_names": [ + "Structure" + ], + "analysis_module": "KratosMultiphysics.StructuralMechanicsApplication", + "analysis_type": "StructuralMechanicsAnalysis", + "analysis_settings": { + "@include_json": "auxiliary_files/system_identification/primal_project_parameters.json" + } + } + } + ], + "responses": [ + { + "name": "damage_response", + "type": "damage_detection_response", + "module": "KratosMultiphysics.SystemIdentificationApplication.responses", + "settings": { + "evaluated_model_part_names": [ + "AdjointStructure" + ], + "adjoint_parameters": { + "@include_json": "auxiliary_files/system_identification_p_norm/adjoint_project_parameters.json" + }, + "test_analysis_list": [ + { + "primal_analysis_name": "Structure_static", + "sensor_measurement_csv_file": "auxiliary_files/damaged_problem/measured_data_ref.csv", + "weight": 1.0 + } + ] + } + } + ], + "controls": [ + { + "name": "material_control", + "type": "material_properties_control", + "module": "KratosMultiphysics.SystemIdentificationApplication.controls", + "settings": { + "model_part_names": [ + { + "primal_model_part_name": "Structure", + "adjoint_model_part_name": "AdjointStructure" + } + ], + "control_variable_name": "YOUNG_MODULUS", + "control_variable_bounds": [ + 0.0, + 30000000000.0 + ], + "filter_settings": { + "filter_type": "explicit_filter", + "filter_radius_settings": { + "filter_radius": 5.0, + "filter_radius_type": "constant" + } + } + } + } + ], + "algorithm_settings": { + "type": "algorithm_steepest_descent", + "settings": { + "echo_level": 0, + "line_search": { + "type": "const_step", + "init_step": 0.01, + "gradient_scaling": "inf_norm" + }, + "conv_settings": { + "type": "max_iter", + "max_iter": 5 + } + }, + "controls": [ + "material_control" + ], + "objective": { + "response_name": "damage_response", + "type": "minimization", + "scaling": 1.0 + } + }, + "processes": { + "kratos_processes": {}, + "optimization_data_processes": { + "output_processes": [ + { + "type": "optimization_problem_ascii_output_process", + "module": "KratosMultiphysics.OptimizationApplication.processes", + "settings": { + "output_file_name": "auxiliary_files/summary_p_norm.csv", + "write_kratos_version": false, + "write_time_stamp": false + } + } + ] + } + } +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/sensor_data.json b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/sensor_data.json new file mode 100644 index 00000000000..3cf0a099f65 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm/sensor_data.json @@ -0,0 +1,109 @@ +{ + "list_of_sensors": [ + { + "type": "displacement_sensor", + "name": "disp_x_99", + "location": [ + 58.54103333333333, + 15.060899999999998, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + }, + { + "type": "displacement_sensor", + "name": "disp_x_108", + "location": [ + 59.34613333333333, + 26.242666666666665, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + }, + { + "type": "displacement_sensor", + "name": "disp_x_89", + "location": [ + 59.49563333333333, + 3.304283333333333, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + }, + { + "type": "displacement_sensor", + "name": "disp_x_422", + "location": [ + 10.403176666666667, + 15.740466666666666, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + }, + { + "type": "displacement_sensor", + "name": "disp_x_142", + "location": [ + 34.94099999999999, + 12.744933333333332, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + }, + { + "type": "displacement_sensor", + "name": "disp_x_486", + "location": [ + 57.3192, + 20.114866666666664, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + }, + { + "type": "displacement_sensor", + "name": "disp_x_94", + "location": [ + 59.23526666666667, + 8.81829, + 0.0 + ], + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "weight": 1.0 + } + ] +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm_summary_ref.csv b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm_summary_ref.csv new file mode 100644 index 00000000000..971304c126c --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm_summary_ref.csv @@ -0,0 +1,20 @@ +# Optimization problem ascii output +# Kratos version: not_given +# Timestamp : not_specified +# ----------------------------------------------- +# --------------- Initial values ---------------- +# algorithm: +# damage_response: +# initial_value: 2.148297976e-04 +# ------------ End of initial values ------------ +# ----------------------------------------------- +# ------------ Start of step values ------------- +# Headers: +# STEP, algorithm:std_obj_value, algorithm:rel_obj[%], algorithm:abs_obj[%], damage_response:value + 0, 2.148297976e-04, 0.000000000e+00, 0.000000000e+00, 2.148297976e-04 + 1, 2.055824507e-04, -4.304499192e+00, -4.304499192e+00, 2.055824507e-04 + 2, 1.966577468e-04, -4.341179814e+00, -8.458812956e+00, 1.966577468e-04 + 3, 1.880459899e-04, -4.379058066e+00, -1.246745469e+01, 1.880459899e-04 + 4, 1.797379508e-04, -4.418088938e+00, -1.633472039e+01, 1.797379508e-04 + 5, 1.717248431e-04, -4.458216926e+00, -2.006470005e+01, 1.717248431e-04 +# End of file \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/test_system_identification.py b/applications/SystemIdentificationApplication/tests/test_system_identification.py index 5e33fc47366..a8a16e95679 100644 --- a/applications/SystemIdentificationApplication/tests/test_system_identification.py +++ b/applications/SystemIdentificationApplication/tests/test_system_identification.py @@ -39,6 +39,25 @@ def test_SystemIdentification(self): }""") CompareTwoFilesCheckProcess(params).Execute() + def test_SystemIdentificationPNorm(self): + model = Kratos.Model() + with open("auxiliary_files/system_identification_p_norm/optimization_parameters.json", "r") as file_input: + params = Kratos.Parameters(file_input.read()) + + analysis = OptimizationAnalysis(model, params) + analysis.Run() + + params = Kratos.Parameters("""{ + "reference_file_name" : "auxiliary_files/system_identification_p_norm_summary_ref.csv", + "output_file_name" : "auxiliary_files/summary_p_norm.csv", + "remove_output_file" : true, + "comparison_type" : "csv_file", + "tolerance" : 1e-5, + "relative_tolerance" : 1e-6, + "dimension" : 3 + }""") + CompareTwoFilesCheckProcess(params).Execute() + if __name__ == '__main__': UnitTest.main()