diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 53773a6f61..9f7a1c913a 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/opm/models/discretization/common/fvbaselinearizer.hh b/opm/models/discretization/common/fvbaselinearizer.hh index 4672091313..d9cf6e0977 100644 --- a/opm/models/discretization/common/fvbaselinearizer.hh +++ b/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; }; diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index f6e7c1d4bc..3a17a705af 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -47,6 +47,7 @@ #include #include #include +#include #include #include @@ -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;igetsize();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;jgetsize();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(A.N()+1); - auto cols = std::make_unique(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;igetsize();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) diff --git a/opm/simulators/linalg/exportSystem.hpp b/opm/simulators/linalg/exportSystem.hpp new file mode 100644 index 0000000000..b42eeb62b7 --- /dev/null +++ b/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 . +*/ + +#ifndef OPM_EXPORT_SYSTEM_HEADER_INCLUDED +#define OPM_EXPORT_SYSTEM_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + + // Forward declaring as they will be called from exportSystem(). + template + void exportSparsity(const IstlMatrix& A, const char *path="."); + template + void exportNonzeros(const IstlMatrix& A, const char *tag="", const char *path="."); + template + void exportVector(const GlobalEqVector& x, const char *tag="", const char *name="export/x"); + + /*! + * \brief Export blocks-sparse linear system. + */ + template + 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 + 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 + 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 + void exportSparsity(const IstlMatrix& A, const char *path) + { + //assemble csr graph + auto rows = std::make_unique(A.N()+1); + auto cols = std::make_unique(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;igetsize();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