[SIAPP]Reverting P norm PR

applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.cpp

@@ -6,8 +6,7 @@
 //
 //  License:         SystemIdentificationApplication/license.txt
 //
-//  Main authors:    Suneth Warnakulasuriya,
-//                   Ihar Antonau
+//  Main authors:    Suneth Warnakulasuriya
 //
 
 // System includes
@@ -90,7 +89,7 @@ struct PartialSensitivity
 
 MeasurementResidualResponseFunction::MeasurementResidualResponseFunction(const double PCoefficient)
     : mPCoefficient(PCoefficient)
-{
+    {
     mResponseGradientList.resize(ParallelUtilities::GetNumThreads());
 }
 
@@ -146,15 +145,19 @@ double MeasurementResidualResponseFunction::CalculateValue(ModelPart& rModelPart
 {
     KRATOS_TRY
 
+    double sum_B_p = 0.0;
     double value = 0.0;
     for (auto& p_sensor : mpSensorsList) {
         const double sensor_value = p_sensor->CalculateValue(rModelPart);
         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);
+        const double current_sensor_error = sensor_value - p_sensor->GetValue(SENSOR_MEASURED_VALUE);
+        p_sensor->SetValue(SENSOR_ERROR, current_sensor_error);
+        sum_B_p += ( std::pow( 0.5 * pow(current_sensor_error, 2) * p_sensor->GetWeight(), mPCoefficient ) );
+        value += p_sensor->GetWeight() * std::pow(sensor_value - p_sensor->GetValue(SENSOR_MEASURED_VALUE), 2) * 0.5;
     }
-    return std::pow(value, 1 / mPCoefficient);
+    mC1 = std::pow( sum_B_p, 1/mPCoefficient - 1 ) / std::pow(2, mPCoefficient - 1);  
+
+    return std::pow(sum_B_p, 1 / mPCoefficient);
 
 
     KRATOS_CATCH("");
@@ -175,26 +178,10 @@ 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...);
-        const double error = std::sqrt( p_sensor->GetValue(SENSOR_ERROR) );
-        noalias(rResponseGradient) += c1 * c2 * error *  local_sensor_response_gradient;
+        noalias(rResponseGradient) += mC1 * (std::pow(p_sensor->GetWeight(), mPCoefficient) * std::pow(p_sensor->GetValue(SENSOR_ERROR), mPCoefficient * 2 - 1 ) ) * local_sensor_response_gradient ;      
     }
-
 
     KRATOS_CATCH("");
 }

applications/SystemIdentificationApplication/custom_sensors/measurement_residual_response_function.h

@@ -6,7 +6,8 @@
 //
 //  License:         SystemIdentificationApplication/license.txt
 //
-//  Main authors:    Suneth Warnakulasuriya
+//  Main authors:    Suneth Warnakulasuriya,
+//                   Ihar Antonau
 //
 
 #pragma once
@@ -157,6 +158,8 @@ class KRATOS_API(DIGITAL_TWIN_APPLICATION) MeasurementResidualResponseFunction :
 
     double mPCoefficient;
 
+    double mC1;
+
     std::vector<Sensor::Pointer> mpSensorsList;
 
     std::vector<Vector> mResponseGradientList;

applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_p_norm_summary_ref.csv

@@ -12,9 +12,9 @@
 # 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
+      1,         2.050244372e-04,     -4.564245970e+00,     -4.564245970e+00,       2.050244372e-04
+      2,         1.955625237e-04,     -4.615017463e+00,     -8.968622685e+00,       1.955625237e-04
+      3,         1.864375960e-04,     -4.665989922e+00,     -1.321613758e+01,       1.864375960e-04
+      4,         1.776430951e-04,     -4.717128485e+00,     -1.730984387e+01,       1.776430951e-04
+      5,         1.691723537e-04,     -4.768404519e+00,     -2.125284501e+01,       1.691723537e-04
 # End of file

applications/SystemIdentificationApplication/tests/auxiliary_files/system_identification_summary_ref.csv

@@ -3,6 +3,7 @@
 # Timestamp     : not_specified
 # -----------------------------------------------
 # --------------- Initial values ----------------
+# 	algorithm:
 # 	damage_response:
 # 		initial_value: 7.924977798e-04
 # ------------ End of initial values ------------
@@ -11,9 +12,9 @@
 # Headers:
 #  STEP, algorithm:std_obj_value, algorithm:rel_obj[%], algorithm:abs_obj[%], damage_response:value
       0,         7.924977798e-04,      0.000000000e+00,      0.000000000e+00,       7.924977798e-04
-      1,         7.578660731e-04,     -4.369943684e+00,     -4.369943684e+00,       7.578660731e-04
-      2,         7.244679939e-04,     -4.406857674e+00,     -8.584224159e+00,       7.244679939e-04
-      3,         6.922659748e-04,     -4.444919496e+00,     -1.264758180e+01,       6.922659748e-04
-      4,         6.612241594e-04,     -4.484087994e+00,     -1.656454110e+01,       6.612241594e-04
-      5,         6.313082950e-04,     -4.524315075e+00,     -2.033942414e+01,       6.313082950e-04
+      1,         7.578660738e-04,     -4.369943603e+00,     -4.369943603e+00,       7.578660738e-04
+      2,         7.244679929e-04,     -4.406857889e+00,     -8.584224288e+00,       7.244679929e-04
+      3,         6.922659728e-04,     -4.444919640e+00,     -1.264758206e+01,       6.922659728e-04
+      4,         6.612241562e-04,     -4.484088173e+00,     -1.656454150e+01,       6.612241562e-04
+      5,         6.313082940e-04,     -4.524314774e+00,     -2.033942428e+01,       6.313082940e-04
 # End of file

applications/SystemIdentificationApplication/tests/test_SystemIdentificationApplication.py

@@ -4,6 +4,7 @@
 import test_adjoint_sensors
 import test_sensor_output_process
 import test_system_identification
+import test_damage_response
 import test_smooth_clamper
 
 def AssembleTestSuites():
@@ -15,6 +16,7 @@ def AssembleTestSuites():
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_sensor_output_process.TestSensorOutputProcess]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_smooth_clamper.TestSmoothClamper]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_system_identification.TestSystemIdentification]))
+    smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_damage_response.TestDamageDetectionResponseBase]))
 
     nightSuite = suites['nightly']
     nightSuite.addTests(smallSuite)

applications/SystemIdentificationApplication/tests/test_damage_response.py

@@ -0,0 +1,79 @@
+from abc import ABC, abstractmethod
+import numpy
+
+import KratosMultiphysics as Kratos
+import KratosMultiphysics.OptimizationApplication as KratosOA
+from KratosMultiphysics.OptimizationApplication.responses.response_routine import ResponseRoutine
+
+import KratosMultiphysics.KratosUnittest as kratos_unittest
+from KratosMultiphysics.OptimizationApplication.optimization_analysis import OptimizationAnalysis
+from KratosMultiphysics.OptimizationApplication.utilities.optimization_problem import OptimizationProblem
+from KratosMultiphysics.SystemIdentificationApplication.responses.damage_detection_response import DamageDetectionResponse
+from KratosMultiphysics.OptimizationApplication.controls.master_control import MasterControl
+
+
+
+class TestDamageDetectionResponseBase(kratos_unittest.TestCase, ABC):
+    def test_damage_response(self):
+        with kratos_unittest.WorkFolderScope(".", __file__):
+            with open("auxiliary_files/system_identification/optimization_parameters.json", "r") as file_input:
+                parameters = Kratos.Parameters(file_input.read())
+
+            model = Kratos.Model()
+            analysis = OptimizationAnalysis(model, parameters)
+
+            analysis.Initialize()
+            analysis.Check()
+            objective: ResponseRoutine = analysis.optimization_problem.GetComponent("damage_response", ResponseRoutine)
+            var = objective.GetRequiredPhysicalGradients()
+            response = analysis.optimization_problem.GetResponse("damage_response")
+            ref_value = response.CalculateValue()
+            self.assertAlmostEqual(ref_value, 0.0007924977797682586, 12)
+            sensitivity = analysis.optimization_problem.GetComponent("master_control", MasterControl).GetEmptyField()
+            response.CalculateGradient(var)
+            gradients = var[Kratos.YOUNG_MODULUS].Evaluate()
+
+            model_part = response.GetInfluencingModelPart()
+
+            for index, element in enumerate(model_part.Elements):
+                element.Properties[Kratos.YOUNG_MODULUS] += 1e3
+                new_value = response.CalculateValue()
+                sensitivity = (new_value - ref_value) / 1e3
+                self.assertAlmostEqual(gradients[index], sensitivity, 12)
+                element.Properties[Kratos.YOUNG_MODULUS] -= 1e3
+                if index == 9:
+                    break
+
+    def test_damage_response_p_norm(self):
+        with kratos_unittest.WorkFolderScope(".", __file__):
+            with open("auxiliary_files/system_identification_p_norm/optimization_parameters.json", "r") as file_input:
+                parameters = Kratos.Parameters(file_input.read())
+
+            model = Kratos.Model()
+            analysis = OptimizationAnalysis(model, parameters)
+
+            analysis.Initialize()
+            analysis.Check()
+            objective: ResponseRoutine = analysis.optimization_problem.GetComponent("damage_response", ResponseRoutine)
+            var = objective.GetRequiredPhysicalGradients()
+            print(var)
+            response = analysis.optimization_problem.GetResponse("damage_response")
+            ref_value = response.CalculateValue()
+            self.assertAlmostEqual(ref_value, 0.00021482979760591695, 12)
+            sensitivity = analysis.optimization_problem.GetComponent("master_control", MasterControl).GetEmptyField()
+            response.CalculateGradient(var)
+            gradients = var[Kratos.YOUNG_MODULUS].Evaluate()
+
+            model_part = response.GetInfluencingModelPart()
+
+            for index, element in enumerate(model_part.Elements):
+                element.Properties[Kratos.YOUNG_MODULUS] += 1e3
+                new_value = response.CalculateValue()
+                sensitivity = (new_value - ref_value) / 1e3
+                self.assertAlmostEqual(gradients[index], sensitivity, 12)
+                element.Properties[Kratos.YOUNG_MODULUS] -= 1e3
+                if index == 9:
+                    break
+
+if __name__ == "__main__":
+    kratos_unittest.main()