Merge pull request #6677 from atgeirr/fix-sprintf

Move linear system export and make it safer.

CMakeLists_files.cmake

@@ -1041,6 +1041,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/DILU.hpp
   opm/simulators/linalg/domesticoverlapfrombcrsmatrix.hh
   opm/simulators/linalg/elementborderlistfromgrid.hh
+  opm/simulators/linalg/exportSystem.hpp
   opm/simulators/linalg/extractMatrix.hpp
   opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp
   opm/simulators/linalg/ExtraSmoothers.hpp

opm/models/discretization/common/fvbaselinearizer.hh

@@ -287,52 +287,6 @@ public:
     GlobalEqVector& residual()
     { return residual_; }
 
-    void printVector(GlobalEqVector&, const char* /*name*/="x")
-    {
-        return;
-    }
-
-    void printResidual(const char* /*name*/="r")
-    {
-        printVector(residual_);
-    }
-
-    void printSparsity(const char* /*name*/="s")
-    {
-        return;
-    }
-
-    void printNonzeros(const char* /*name*/="d")
-    {
-        return;
-    }
-
-    void printJacobian()
-    {
-        return;
-    }
-
-    void exportSystem(int /*idx*/, char* /*tag*/, const char* /*path*/="export")
-    {
-        return;
-    }
-
-    void exportVector(GlobalEqVector& x, const char* /*tag*/="", const char* /*name*/="export/x")
-    {
-        printf("n = %lu\n",x.dim());
-    }
-
-    void exportSparsity(const char* /*path*/=".")
-    {
-        return;
-    }
-
-    void exportNonzeros(const char* /*tag*/="", const char* /*path*/=".")
-    {
-        return;
-    }
-
-
     void setLinearizationType(LinearizationType linearizationType){
         linearizationType_ = linearizationType;
     };

opm/models/discretization/common/tpfalinearizer.hh

@@ -47,6 +47,7 @@
 #include <opm/models/discretization/common/baseauxiliarymodule.hh>
 #include <opm/models/discretization/common/fvbaseproperties.hh>
 #include <opm/models/discretization/common/linearizationtype.hh>
+#include <opm/simulators/linalg/exportSystem.hpp>
 
 #include <cassert>
 #include <cstddef>
@@ -419,198 +420,21 @@ public:
     { return residual_; }
 
     /*!
-     * \brief Print first 16 block elements of a vector.
+     * \brief Export block sparse linear system.
      */
-    void printVector(GlobalEqVector &x, const char *name="x")
+    void exportSystem(const int idx, std::string& tag, const char *path="export")
     {
-        int count = 1;
-        fmt::print("{} =\n[\n",name);
-        for (auto block = x.begin(); block != x.end(); block++)
-        {
-            for (auto i=block->begin(); i!=block->end(); i++)
-            {
-                fmt::print(" {:+.4e}",*i);
-            }
-            fmt::print("\n");
-            count++;
-            if(count>16) break;
-        }
-        fmt::print("]\n");
-    }
-
-    /*!
-     * \brief Print first 16 block elements of residual.
-     */
-    void printResidual(const char *name="r")
-    {
-        printVector(residual_, name);
-    }
-
-    /*!
-     * \brief Print sparsity pattern of the first 16 rows
-     *        of the jacobian block matrix
-     */
-    void printSparsity(const char *name="s")
-    {
-        auto& A = jacobian_->istlMatrix();
-
-        fmt::print("nrows = {:d}\n",A.N());
-        fmt::print("ncols = {:d}\n",A.M());
-        fmt::print("nnz   = {:d}\n",A.nonzeroes());
-
-        fmt::print("{} =\n[\n",name);
-        int count=1;
-        int offset=0;
-        for(auto row=A.begin(); row!=A.end(); row++)
-        {
-            fmt::print("{:4d}: ",offset);
-            for(unsigned int i=0;i<row->getsize();i++)
-            {
-                fmt::print(" {:4d}",row->getindexptr()[i]);
-            }
-            fmt::print("\n");
-            offset+=row->getsize();
-            count++;
-            if(count>16) break;
-        }
-        fmt::print("]\n");
-    }
-
-    /*!
-     * \brief Print block elements of the first 6 rows of the
-     * j      acobian block matrix
-     */
-    void printNonzeros(const char *name="d")
-    {
-        auto& A = jacobian_->istlMatrix();
-        fmt::print("{} =\n[\n",name);
-        int count=1;
-        for(auto row=A.begin();row!=A.end();row++)
-        {
-            for(unsigned int j=0;j<row->getsize();j++)
-            {
-                fmt::print("|");
-                auto mat = row->getptr()[j];
-                for(auto vec=mat.begin();vec!=mat.end();vec++)
-                {
-                    for(auto k=vec->begin();k!=vec->end();k++)
-                    {
-                        fmt::print(" {:+.4e}",*k);
-                    }
-                    fmt::print(" |");
-                }
-                fmt::print("\n");
-            }
-            count++;
-            if(count>6) break;
-            fmt::print("\n");
-        }
-        fmt::print("]\n");
-    }
-
-    /*!
-     * \brief Print sparsity pattern and nonzeros of jacobian block matrix
-     */
-    void printJacobian()
-    {
-        printSparsity();
-        printNonzeros();
-    }
-
-    /*!
-     * \brief Export blocks-sparse linear system.
-     */
-    void exportSystem(int idx, char *tag, const char *path="export")
-    {
-        // export sparsity only once
-        if(exportIndex_==-1) exportSparsity(path);
+        const bool export_sparsity = exportIndex_ == -1;
 
         // increment indices and generate tag
-        exportCount_ = exportIndex_==idx ? ++exportCount_ : 0;
+        exportCount_ = exportIndex_ == idx ? ++exportCount_ : 0;
         exportIndex_ = idx;
-        sprintf(tag,"_%03d_%02d",exportIndex_, exportCount_);
+        tag = fmt::format("_{:03d}_{:02d}", exportIndex_, exportCount_);
 
         fmt::print("index = {:d}\n", exportIndex_);
         fmt::print("count = {:d}\n", exportCount_);
 
-        // export matrix
-        exportNonzeros(tag,path);
-
-        // export residual
-        char name[256];
-        sprintf(name,"%s/r",path);
-        exportVector(residual_,tag,name);
-    }
-
-    /*!
-     * \brief Export block vector.
-     */
-    void exportVector(GlobalEqVector &x, const char *tag="", const char *name="export/x")
-    {
-        // assume double precision and contiguous data
-        const double *data = &x[0][0];
-
-        char filename[512];
-        sprintf(filename,"%s%s.f64",name,tag);
-        FILE *out =fopen(filename,"w");
-        fwrite(data, sizeof(double), x.dim(),out);
-        fclose(out);
-    }
-
-    /*!
-     * \brief Export nonzero blocks of jacobian block-sparse matrix
-     */
-    void exportNonzeros(const char *tag="", const char *path=".")
-    {
-        auto& A = jacobian_->istlMatrix();
-
-        // assume double precision and contiguous data
-        const double *data = &A[0][0][0][0];
-        size_t dim = A[0][0].N()*A[0][0].M()*A.nonzeroes();
-
-        char filename[256];
-        sprintf(filename,"%s/data%s.f64",path,tag);
-        FILE *out =fopen(filename,"w");
-        fwrite(data, sizeof(double), dim,out);
-        fclose(out);
-    }
-
-    /*!
-     * \brief Export sparsity pattern of jacobian block-sparse matrix
-     */
-    void exportSparsity(const char *path=".")
-    {
-        //assemble csr graph
-        auto& A = jacobian_->istlMatrix();
-        auto rows = std::make_unique<int[]>(A.N()+1);
-        auto cols = std::make_unique<int[]>(A.nonzeroes());
-
-        int irow=0;
-        int icol=0;
-        rows[0]=0;
-        for(auto row=A.begin(); row!=A.end(); row++)
-        {
-            for(unsigned int i=0;i<row->getsize();i++)
-            {
-                cols[icol++]=row->getindexptr()[i];
-            }
-            rows[irow+1]= rows[irow]+row->getsize();
-            irow++;
-        }
-
-        //export arrays
-        FILE *out;
-        char filename[256];
-
-        sprintf(filename,"%s/rows.i32",path);
-        out=fopen(filename,"w");
-        fwrite(rows, sizeof(int), A.N()+1,out);
-        fclose(out);
-
-        sprintf(filename,"%s/cols.i32",path);
-        out=fopen(filename,"w");
-        fwrite(cols, sizeof(int), A.nonzeroes(),out);
-        fclose(out);
+        Opm::exportSystem(jacobian_->istlMatrix(), residual_, export_sparsity, tag.c_str(), path);
     }
 
     void setLinearizationType(LinearizationType linearizationType)

opm/simulators/linalg/exportSystem.hpp

@@ -0,0 +1,139 @@
+/*
+  Copyright 2025 Equinor ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_EXPORT_SYSTEM_HEADER_INCLUDED
+#define OPM_EXPORT_SYSTEM_HEADER_INCLUDED
+
+#include <cstdio>
+#include <memory>
+
+namespace Opm
+{
+
+    // Forward declaring as they will be called from exportSystem().
+    template <class IstlMatrix>
+    void exportSparsity(const IstlMatrix& A, const char *path=".");
+    template <class IstlMatrix>
+    void exportNonzeros(const IstlMatrix& A, const char *tag="", const char *path=".");
+    template <class GlobalEqVector>
+    void exportVector(const GlobalEqVector& x, const char *tag="", const char *name="export/x");
+
+    /*!
+     * \brief Export blocks-sparse linear system.
+     */
+    template <class IstlMatrix, class GlobalEqVector>
+    void
+    exportSystem(const IstlMatrix& jacobian,
+                 const GlobalEqVector& residual,
+                 const bool export_sparsity,
+                 const char* tag,
+                 const char* path = "export")
+    {
+        // export sparsity only if requested
+        if (export_sparsity) {
+            exportSparsity(jacobian,path);
+        }
+
+        // export matrix
+        exportNonzeros(jacobian,tag,path);
+
+        // export residual
+        constexpr size_t bufsize = 256;
+        char name[bufsize];
+        std::snprintf(name,bufsize,"%s/r",path);
+        exportVector(residual,tag,name);
+    }
+
+    /*!
+     * \brief Export block vector.
+     */
+    template <class GlobalEqVector>
+    void exportVector(const GlobalEqVector& x, const char *tag, const char *name)
+    {
+        // assume double precision and contiguous data
+        const double *data = &x[0][0];
+
+        constexpr size_t bufsize = 512;
+        char filename[bufsize];
+        std::snprintf(filename,bufsize,"%s%s.f64",name,tag);
+        FILE *out =fopen(filename,"w");
+        std::fwrite(data, sizeof(double), x.dim(),out);
+        std::fclose(out);
+    }
+
+    /*!
+     * \brief Export nonzero blocks of jacobian block-sparse matrix
+     */
+    template <class IstlMatrix>
+    void exportNonzeros(const IstlMatrix& A, const char *tag, const char *path)
+    {
+        // assume double precision and contiguous data
+        const double *data = &A[0][0][0][0];
+        size_t dim = A[0][0].N()*A[0][0].M()*A.nonzeroes();
+
+        constexpr size_t bufsize = 256;
+        char filename[bufsize];
+        std::snprintf(filename,bufsize,"%s/data%s.f64",path,tag);
+        FILE *out =fopen(filename,"w");
+        std::fwrite(data, sizeof(double), dim,out);
+        std::fclose(out);
+    }
+
+    /*!
+     * \brief Export sparsity pattern of jacobian block-sparse matrix
+     */
+    template <class IstlMatrix>
+    void exportSparsity(const IstlMatrix& A, const char *path)
+    {
+        //assemble csr graph
+        auto rows = std::make_unique<int[]>(A.N()+1);
+        auto cols = std::make_unique<int[]>(A.nonzeroes());
+
+        int irow=0;
+        int icol=0;
+        rows[0]=0;
+        for(auto row=A.begin(); row!=A.end(); row++)
+        {
+            for(unsigned int i=0;i<row->getsize();i++)
+            {
+                cols[icol++]=row->getindexptr()[i];
+            }
+            rows[irow+1]= rows[irow]+row->getsize();
+            irow++;
+        }
+
+        //export arrays
+        FILE *out;
+        constexpr size_t bufsize = 256;
+        char filename[bufsize];
+
+        std::snprintf(filename,bufsize,"%s/rows.i32",path);
+        out=std::fopen(filename,"w");
+        std::fwrite(rows.get(), sizeof(int), A.N()+1,out);
+        std::fclose(out);
+
+        std::snprintf(filename,bufsize,"%s/cols.i32",path);
+        out=std::fopen(filename,"w");
+        std::fwrite(cols.get(), sizeof(int), A.nonzeroes(),out);
+        std::fclose(out);
+    }
+
+} // namespace Opm
+
+#endif // OPM_EXPORT_SYSTEM_HEADER_INCLUDED