SAR: refactor and improve SAR calculations

* new range limit with auto range option to reduce required computation
* averaging done only once for all frequencies (massive speedup)
* split read and write hdf5 into separate functions
* use the (now separate) hdf5 write in the openEMS field proc to reduce
  redundancy

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

Author: thliebig

Common/processfields_sar.cpp

@@ -272,46 +272,36 @@ void ProcessFieldsSAR::DumpFDData()
 		SAR_Calc.SetCellDensities(&cell_density);
 		SAR_Calc.SetCellWidth(cellWidth);
 		SAR_Calc.SetCellVolumes(&cell_volume);
-		SAR_Calc.SetCellCondictivity(&cell_kappa); // cell_kappa will be NULL if m_UseCellKappa is false
-		double mass = SAR_Calc.CalcTotalMass();
 
 		for (size_t n = 0; n<m_FD_Samples.size(); ++n)
 		{
-			SAR_Calc.SetEField(m_E_FD_Fields.at(n));
 			if (!m_UseCellKappa)
-				SAR_Calc.SetJField(m_J_FD_Fields.at(n));
-			power = SAR_Calc.CalcSARPower();
-			SAR_Calc.SetCellVolumes(&cell_volume);
-			power = SAR_Calc.CalcSARPower();
-			SAR_Calc.CalcSAR(SAR);
+				SAR_Calc.AddEFieldAndJField(m_FD_Samples.at(n), m_E_FD_Fields.at(n), m_J_FD_Fields.at(n));
+			else
+				SAR_Calc.AddEFieldAndCondictivity(m_FD_Samples.at(n), m_E_FD_Fields.at(n), &cell_kappa);
+		}
+		SAR_Calc.CalcSAR();
 
-			if (m_fileType==VTK_FILETYPE)
+		if (m_fileType==VTK_FILETYPE)
+		{
+			for (size_t n = 0; n<m_FD_Samples.size(); ++n)
 			{
 				stringstream ss;
 				ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n);
-
 				m_Vtk_Dump_File->SetFilename(ss.str());
 				m_Vtk_Dump_File->ClearAllFields();
-				m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR);
+				m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType), *SAR_Calc.GetSAR(n));
 				if (m_Vtk_Dump_File->Write()==false)
 					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
 			}
-			else if (m_fileType==HDF5_FILETYPE)
-			{
-				stringstream ss;
-				ss << "f" << n;
-				if (m_HDF5_Dump_File->WriteScalarField<float>(ss.str(), SAR, g_settings.GetLegacyHFD5Dumps())==false)
-					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
-				if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",(float)m_FD_Samples.at(n))==false)
-					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
-				if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"mass",mass)==false)
-					cerr << "ProcessFieldsSAR::DumpFDData: can't dump mass to file...! " << endl;
-				if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",power)==false)
-					cerr << "ProcessFieldsSAR::DumpFDData: can't dump power to file...! " << endl;
-			}
-			else
-				cerr << "ProcessFieldsSAR::DumpFDData: unknown File-Type" << endl;
 		}
+		else if (m_fileType==HDF5_FILETYPE)
+		{
+			if (SAR_Calc.WriteToHDF5(*m_HDF5_Dump_File, g_settings.GetLegacyHFD5Dumps())==false)
+				cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+		}
+		else
+			cerr << "ProcessFieldsSAR::DumpFDData: unknown File-Type" << endl;
 	}
 	for (int n=0;n<3;++n)
 		delete[] cellWidth[n];

python/openEMS/sar_calculation.pxd

@@ -27,7 +27,9 @@ cdef extern from "openEMS/sar_calculation.h":
         void EnableProgress(bool enable)
         void SetAveragingMass(float mass)
         bool SetAveragingMethod(string method, bool silent)
-        bool CalcFromHDF5(string h5_fn, string out_name, bool export_cube_stats) nogil
+        void EnableCubeStats()
+        void EnableAutoRange(double dBmax)
+        bool CalcFromHDF5(string h5_fn, string out_name) nogil
 
 
 cdef class SAR_Calculation:

python/openEMS/sar_calculation.pyx

@@ -31,6 +31,11 @@ cdef class SAR_Calculation:
                 self.SetDebugLevel(int(v))
             elif k=='debug':
                 self.SetDebugLevel(int(v))
+            elif k=='autoRange':
+                self.EnableAutoRange(float(v))
+            elif k=='EnableCubeStats':
+                if v:
+                    self.EnableCubeStats()
             else:
                 raise Exception('Unknown keyword argument: "{}"'.format(k))
 
@@ -49,13 +54,20 @@ cdef class SAR_Calculation:
          #cdef string c_method = method.encode('UTF-8')
          return self.thisptr.SetAveragingMethod(method.encode('UTF-8'), silent)
 
+    def EnableAutoRange(self, dBmax):
+        self.thisptr.EnableAutoRange(float(dBmax))
+
+    def EnableCubeStats(self):
+        self.thisptr.EnableCubeStats()
+
     def CalcFromHDF5(self, h5_fn, out_name, export_cube_stats=False):
         if not os.path.exists(h5_fn):
             raise Exception('File "{}" does not exist'.format(h5_fn))
         cdef string in_fn = h5_fn.encode('UTF-8')
         cdef string out_fn = out_name.encode('UTF-8')
-        cdef bool ecs = export_cube_stats
+        if export_cube_stats:
+            self.EnableCubeStats()
         with nogil:
-            ok = self.thisptr.CalcFromHDF5(in_fn, out_fn, ecs)
+            ok = self.thisptr.CalcFromHDF5(in_fn, out_fn)
         return ok
 

tools/CMakeLists.txt

@@ -35,7 +35,7 @@ ADD_EXECUTABLE( sar_bin sar_calc.cpp)
 SET_TARGET_PROPERTIES(sar_bin PROPERTIES OUTPUT_NAME sar_calc)
 TARGET_LINK_LIBRARIES(sar_bin openEMS)
 
-INSTALL(FILES sar_calculation.h DESTINATION include/openEMS)
+INSTALL(FILES sar_calculation.h hdf5_file_reader.h hdf5_file_writer.h vtk_file_writer.h DESTINATION include/openEMS)
 
 # not ideal, but no way around it?
 INSTALL(FILES arraylib/array_ij.h arraylib/array_ijk.h arraylib/array_nijk.h DESTINATION include/openEMS/arraylib)

tools/hdf5_file_writer.cpp

@@ -420,15 +420,15 @@ bool HDF5_File_Writer::WriteAtrribute(hid_t loc, std::string attr_name, void con
 	hid_t attribute_id = H5Acreate(loc, attr_name.c_str(), mem_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
 	if (attribute_id<0)
 	{
-		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to create the attribute" << endl;
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to create the attribute: " << attr_name << endl;
 		H5Sclose(dataspace_id);
 		return false;
 	}
 
 	/* Write the attribute data. */
 	if (H5Awrite(attribute_id, mem_type, value)<0)
 	{
-		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to write the attribute" << endl;
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to write the attribute" << attr_name << endl;
 		H5Aclose(attribute_id);
 		H5Sclose(dataspace_id);
 		return false;
@@ -468,6 +468,26 @@ bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name
 	return ok;
 }
 
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, vector<int> values)
+{
+	int* val = new int[values.size()];
+	for (size_t n=0;n<values.size();++n)
+		val[n]=values.at(n);
+	bool ok =  HDF5_File_Writer::WriteAtrribute(locName, attr_name, val, values.size(), H5T_NATIVE_INT);
+	delete[] val; val=NULL;
+	return ok;
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, vector<unsigned int> values)
+{
+	unsigned int* val = new unsigned int[values.size()];
+	for (size_t n=0;n<values.size();++n)
+		val[n]=values.at(n);
+	bool ok =  HDF5_File_Writer::WriteAtrribute(locName, attr_name, val, values.size(), H5T_NATIVE_UINT);
+	delete[] val; val=NULL;
+	return ok;
+}
+
 bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, float value)
 {
 	return HDF5_File_Writer::WriteAtrribute(locName, attr_name, &value, 0, H5T_NATIVE_FLOAT);

tools/hdf5_file_writer.h

@@ -53,6 +53,8 @@ class HDF5_File_Writer
 	bool WriteAtrribute(std::string locName, std::string attr_name, double const* value, hsize_t size);
 	bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<float> values);
 	bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<double> values);
+	bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<int> values);
+	bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<unsigned int> values);
 	bool WriteAtrribute(std::string locName, std::string attr_name, float value);
 	bool WriteAtrribute(std::string locName, std::string attr_name, double value);
 	bool WriteAtrribute(std::string locName, std::string attr_name, int value);

tools/sar_calc.cpp

@@ -27,14 +27,15 @@ int main(int argc, char *argv[])
 {
 	cout << " ---------------------------------------------------------------------- " << endl;
 	cout << " | SAR calculation for openEMS "                                          << endl;
-	cout << " | (C) 2012-2025 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license"   << endl;
+	cout << " | (C) 2012-2026 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license"   << endl;
 	cout << " ---------------------------------------------------------------------- " << endl;
 
 	options_description desc("Options");
 	variables_map vm;
 	try {
 		string ifile, ofile, method;
 		double m_avg = 0;
+		double auto_range=0;
 		bool debug = false;
 		bool export_cube_stats = false;
 		bool legacyHDF5 = false;
@@ -46,6 +47,7 @@ int main(int argc, char *argv[])
 		("output,o", value(&ofile)->required(), "pathname for output hdf5 file")
 		("method",   value(&method)->default_value("SIMPLE"), "set SAR method: IEEE_C95_3, IEEE_62704, SIMPLE")
 		("mass,m"  , value(&m_avg), "averaging mass in g")
+		("autorange,a"  , value(&auto_range), "autorange, value limit in dB from max. (>0)")
 		("verbose,v", bool_switch(&debug), "verbose")
 		("progress,p", bool_switch(&progress), "show progress")
 		("export_cube_stats,e", bool_switch(&export_cube_stats), "Export Cube Statistics")
@@ -64,11 +66,13 @@ int main(int argc, char *argv[])
 		sar_calc.SetDebugLevel(int(debug));
 		sar_calc.EnableProgress(progress);
 		sar_calc.SetAveragingMass(m_avg/1000);
+		sar_calc.EnableAutoRange(auto_range);
+		if (export_cube_stats)
+			sar_calc.EnableCubeStats();
 		//sar_calc.SetAveragingMethod(SAR_Calculation::IEEE_62704);
 		if (!sar_calc.SetAveragingMethod(method, !debug))
 			return -1;
-		sar_calc.CalcFromHDF5(ifile, ofile, export_cube_stats, legacyHDF5);
-		return 0;
+		return sar_calc.CalcFromHDF5(ifile, ofile, legacyHDF5);
 	}
 	catch(exception& e) {
 		cerr << e.what() << "\n";