Merge branch 'CurveFromFile' into 'master'

Read curve definitions from binary file

See merge request ogs/ogs!5127

Applications/FileIO/GocadIO/GocadSGridReader.cpp

@@ -504,8 +504,8 @@ void GocadSGridReader::readElementPropertiesBinary()
             WARN("Empty filename for property {:s}.", property._property_name);
             continue;
         }
-        std::vector<float> float_properties =
-            BaseLib::readBinaryArray<float>(fname, _index_calculator._n_cells);
+        std::vector<float> float_properties = BaseLib::readBinaryVector<float>(
+            fname, 0, _index_calculator._n_cells);
         DBUG(
             "GocadSGridReader::readElementPropertiesBinary(): Read {:d} float "
             "properties from binary file.",

BaseLib/FileTools.cpp

@@ -17,6 +17,9 @@
 #include <spdlog/fmt/bundled/core.h>
 
 #include <boost/algorithm/string/predicate.hpp>
+#include <boost/endian/conversion.hpp>
+#include <boost/interprocess/file_mapping.hpp>
+#include <boost/interprocess/mapped_region.hpp>
 #include <filesystem>
 #include <fstream>
 #include <typeindex>
@@ -277,6 +280,21 @@ bool createOutputDirectory(std::string const& dir)
     return true;
 }
 
+std::vector<double> readDoublesFromBinaryFile(const std::string& filename)
+{
+    auto prj_dir = BaseLib::getProjectDirectory();
+    std::string path_to_file = BaseLib::joinPaths(prj_dir, filename);
+    std::string file_extension = BaseLib::getFileExtension(filename);
+    if (file_extension != ".bin")
+    {
+        OGS_FATAL(
+            "Currently only binary files with extension '.bin' supported. The "
+            "specified file has extension {:s}.",
+            file_extension)
+    }
+    return BaseLib::readBinaryVector<double>(path_to_file);
+}
+
 template <typename T>
 T readBinaryValue(std::istream& in)
 {
@@ -290,47 +308,62 @@ template float readBinaryValue<float>(std::istream&);
 template double readBinaryValue<double>(std::istream&);
 
 template <typename T>
-std::vector<T> readBinaryArray(std::string const& filename, std::size_t const n)
+std::vector<T> readBinaryVector(std::string const& filename,
+                                std::size_t const start_element,
+                                std::size_t const num_elements)
 {
-    std::ifstream in(filename.c_str());
-    if (!in)
+    if (!IsFileExisting(filename))
     {
-        ERR("readBinaryArray(): Error while reading from file '{:s}'.",
-            filename);
-        ERR("Could not open file '{:s}' for input.", filename);
-        in.close();
-        return std::vector<T>();
+        OGS_FATAL("File {:s} not found", filename);
     }
 
-    std::vector<T> result;
-    result.reserve(n);
+    // Determine file size
+    std::uintmax_t file_size = std::filesystem::file_size(filename);
+    std::size_t total_elements = file_size / sizeof(T);
 
-    for (std::size_t p = 0; in && !in.eof() && p < n; ++p)
+    if (start_element >= total_elements)
     {
-        result.push_back(BaseLib::readBinaryValue<T>(in));
+        OGS_FATAL("Start element is beyond file size");
     }
 
-    if (result.size() == n)
-    {
-        return result;
-    }
+    // Calculate the number of elements to read
+    std::size_t const elements_to_read =
+        std::min(num_elements, total_elements - start_element);
+
+    // Calculate offset and size to map
+    std::size_t const offset = start_element * sizeof(T);
+    std::size_t const size_to_map = elements_to_read * sizeof(T);
+
+    // Create a file mapping
+    boost::interprocess::file_mapping file(filename.c_str(),
+                                           boost::interprocess::read_only);
 
-    ERR("readBinaryArray(): Error while reading from file '{:s}'.", filename);
-    ERR("Read different number of values. Expected {:d}, got {:d}.",
-        n,
-        result.size());
+    // Map the specified region
+    boost::interprocess::mapped_region region(
+        file, boost::interprocess::read_only, offset, size_to_map);
 
-    if (!in.eof())
+    // Get the address of the mapped region
+    auto* addr = region.get_address();
+
+    // Create vector and copy data
+    std::vector<T> result(elements_to_read);
+    std::memcpy(result.data(), addr, size_to_map);
+
+    if constexpr (std::endian::native != std::endian::little)
     {
-        ERR("EOF reached.\n");
+        boost::endian::endian_reverse_inplace(result);
     }
 
-    return std::vector<T>();
+    return result;
 }
 
 // explicit template instantiation
-template std::vector<float> readBinaryArray<float>(std::string const&,
-                                                   std::size_t const);
+template std::vector<float> readBinaryVector<float>(std::string const&,
+                                                    std::size_t const,
+                                                    std::size_t const);
+template std::vector<double> readBinaryVector<double>(std::string const&,
+                                                      std::size_t const,
+                                                      std::size_t const);
 
 template <typename T>
 void writeValueBinary(std::ostream& out, T const& val)

BaseLib/FileTools.h

@@ -16,6 +16,7 @@
 
 #include <cstddef>
 #include <iosfwd>
+#include <limits>
 #include <string>
 #include <tuple>
 #include <vector>
@@ -52,6 +53,8 @@ std::string constructFormattedFileName(std::string const& format_specification,
                                        double const t,
                                        int const iteration);
 
+std::vector<double> readDoublesFromBinaryFile(const std::string& filename);
+
 /**
  * \brief write value as binary into the given output stream
  *
@@ -89,8 +92,9 @@ extern template float readBinaryValue<float>(std::istream&);
 extern template double readBinaryValue<double>(std::istream&);
 
 template <typename T>
-std::vector<T> readBinaryArray(std::string const& filename,
-                               std::size_t const n);
+std::vector<T> readBinaryVector(
+    std::string const& filename, std::size_t const start_element = 0,
+    std::size_t const num_elements = std::numeric_limits<std::size_t>::max());
 
 /**
  * Extracts basename from given pathname with extension.

Documentation/ProjectFile/curve/t_read_from_file.md

@@ -0,0 +1 @@
+The default is false. If this option is true, the coordinates and values of the curve are read from a binary file. The precision of the values is double and the byte ordering must be little endian. This option is recommended for large curve definitions, which may fail to directly read from the project file.

Documentation/ProjectFile/prj/time_loop/output/t_fixed_output_times_from_file.md

@@ -0,0 +1 @@
+Read **\<fixed_output_times\>** from binary file (double precision, little endian).

MathLib/Curve/CreatePiecewiseLinearCurve.cpp

@@ -12,21 +12,50 @@
 
 #include "CreatePiecewiseLinearCurve.h"
 
+#include <boost/endian/conversion.hpp>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
+#include "BaseLib/FileTools.h"
+#include "BaseLib/StringTools.h"
 
 namespace MathLib
 {
-
 PiecewiseLinearCurveConfig parsePiecewiseLinearCurveConfig(
     BaseLib::ConfigTree const& config)
 {
-    auto x =
-        //! \ogs_file_param{curve__coords}
-        config.getConfigParameter<std::vector<double>>("coords");
-    auto y =
-        //! \ogs_file_param{curve__values}
-        config.getConfigParameter<std::vector<double>>("values");
+    const bool read_from_file =  //! \ogs_file_param{curve__read_from_file}
+        config.getConfigParameter<bool>("read_from_file", false);
+
+    std::vector<double> x;
+    std::vector<double> y;
+
+    if (read_from_file == true)
+    {
+        auto const coords_file_name =
+            //! \ogs_file_param{curve__coords}
+            config.getConfigParameter<std::string>("coords");
+
+        auto const values_file_name =
+            //! \ogs_file_param{curve__values}
+            config.getConfigParameter<std::string>("values");
+
+        x = BaseLib::readDoublesFromBinaryFile(coords_file_name);
+
+        y = BaseLib::readDoublesFromBinaryFile(values_file_name);
+    }
+    else
+    {
+        x =
+            //! \ogs_file_param{curve__coords}
+            config.getConfigParameter<std::vector<double>>("coords");
+        y =
+            //! \ogs_file_param{curve__values}
+            config.getConfigParameter<std::vector<double>>("values");
+    }
 
     if (x.empty() || y.empty())
     {

MathLib/Curve/CreatePiecewiseLinearCurve.h

@@ -28,6 +28,8 @@ struct PiecewiseLinearCurveConfig
     std::vector<double> ys;
 };
 
+std::vector<double> readDoublesFromBinaryFile(const std::string& filename);
+
 PiecewiseLinearCurveConfig parsePiecewiseLinearCurveConfig(
     BaseLib::ConfigTree const& config);
 

ProcessLib/HeatTransportBHE/Tests.cmake

@@ -96,6 +96,20 @@ AddTest(
     beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_newton_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-13
 )
 
+AddTest(
+    NAME HeatTransportBHE_1U_3D_beier_sandbox_binary_curve
+    PATH Parabolic/T/3D_Beier_sandbox
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS beier_sandbox_binary_curve.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    RUNTIME 20
+    DIFF_DATA
+    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_binary_curve_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 1e-14
+    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_binary_curve_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-12
+)
+
 AddTest(
     NAME HeatTransportBHE_1U_3D_beier_sandbox_algebraicBC
     PATH Parabolic/T/3D_Beier_sandbox

ProcessLib/Output/CreateOutputConfig.cpp

@@ -14,6 +14,7 @@
 
 #include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTree.h"
+#include "BaseLib/FileTools.h"  // required for reading output_times from binary file
 #include "MaterialLib/Utils/MediaCreation.h"  // required for splitMaterialIDString
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Utils/createMaterialIDsBasedSubMesh.h"
@@ -146,6 +147,18 @@ OutputConfig createOutputConfig(
         //! \ogs_file_param{prj__time_loop__output__fixed_output_times}
         config.getConfigParameter<std::vector<double>>("fixed_output_times",
                                                        {});
+    if (output_config.fixed_output_times.empty())
+    {
+        //! \ogs_file_param{prj__time_loop__output__fixed_output_times_from_file}
+        std::string filename = config.getConfigParameter<std::string>(
+            "fixed_output_times_from_file", "no_file");
+
+        if (filename != "no_file")
+        {
+            output_config.fixed_output_times =
+                BaseLib::readDoublesFromBinaryFile(filename);
+        }
+    }
     // Remove possible duplicated elements and sort.
     BaseLib::makeVectorUnique(output_config.fixed_output_times);
 

Tests/Data/Parabolic/T/3D_Beier_sandbox/beier_sandbox_binary_curve.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="beier_sandbox.prj">
+    <add sel="/*/curves/curve">
+        <read_from_file>true</read_from_file>
+    </add>
+    <replace sel="/*/time_loop/output/prefix/text()">beier_sandbox_binary_curve</replace>
+    <replace sel="/*/curves/curve/coords/text()">inflow_temperature_coords.bin</replace>
+    <replace sel="/*/curves/curve/values/text()">inflow_temperature_values.bin</replace>
+</OpenGeoSysProjectDiff>

Tests/Data/Parabolic/T/3D_Beier_sandbox/read_values_from_prj_curve.py

@@ -0,0 +1,33 @@
+import inspect
+import os
+import xml.etree.ElementTree as ET
+
+import numpy as np
+
+file = inspect.getfile(inspect.currentframe())
+prj_path = os.path.split(os.path.realpath(file))[0]
+
+# read prj-data
+tree = ET.parse(prj_path + "/beier_sandbox.prj")
+root = tree.getroot()
+
+curve_elements = root.findall(".//curve")
+for curve in curve_elements:
+    curve_name = curve.find("name").text
+
+    # read coordinates
+    coords_element = curve.find("coords")
+    coords_text = coords_element.text  # replace("\n", "    ")
+    coords_array = np.fromstring(coords_text, dtype=float, sep="    ")
+    # Store array in binary format
+    coords_array.astype("<f8").tofile(
+        prj_path + f"/{curve_name}_coords.bin"
+    )  # < - little Endian | f - float | 8 - Bytes equal to 64bit(double)
+
+    # read values
+    values_element = curve.find("values")
+    values_text = values_element.text  # replace("\n", "    ")
+    values_array = np.fromstring(values_text, dtype=float, sep="    ")
+    values_array.astype("<f8").tofile(
+        prj_path + f"/{curve_name}_values.bin"
+    )  # < - little Endian | f - float | 8 - Bytes equal to 64bit(double)

web/content/docs/userguide/blocks/curves.md

@@ -47,7 +47,24 @@ Curves can be defined using the following block:
 </curves>
 ```
 
-and can be called in "Properties" and "Parameters" inside "Expression":
+or recommended for long curves definitions, they can be read from a binary file (double precision values in little endian format) placed in the project file directory:
+
+```xml
+<curves>
+    <curve>
+        <name>CurveName</name>
+        <read_from_file>true</read_from_file>
+        <coords>
+            coords_filename.bin
+        </coords>
+        <values>
+            values_filename.bin
+        </values>
+    </curve>
+</curves>
+```
+
+They can be called in "Properties" and "Parameters" inside "Expression":
 `CurveName(evaluation_value)`
 
 where `evaluation_value` is always exactly one value and refers to values provided inside `<coords> </coords>`.