diff --git a/.github/workflows/build_guide.yml b/.github/workflows/build_guide.yml new file mode 100644 index 00000000..d132cb08 --- /dev/null +++ b/.github/workflows/build_guide.yml @@ -0,0 +1,39 @@ +name: build_guide + +# Generate sphinx docs + +on: + workflow_dispatch: + +jobs: + + build_docs_job: + runs-on: ubuntu-latest + permissions: + contents: write + steps: + - name: Checkout + uses: actions/checkout@v2.3.4 + + - name: Set up Python + uses: actions/setup-python@v2.2.1 + with: + python-version: 3.9 + + - name: Install dependencies + run: | + python -m pip install -U sphinx + python -m pip install sphinx-rtd-theme + python -m pip install breathe + python -m pip install ghp-import + + - name: make the sphinx docs + run: | + make -C docs clean + make -C docs html + + - name: Run ghp-import + run: | + ghp-import -n -p -f ~/local/doc/build/html + + diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml index 42421825..098decaf 100644 --- a/.github/workflows/c-cpp.yml +++ b/.github/workflows/c-cpp.yml @@ -1,123 +1,102 @@ name: LBPM CI on: + workflow_dispatch: push: - branches: [ master ] + branches: [master, dev] pull_request: - branches: [ master ] + branches: [master, dev] jobs: build-and-test: - runs-on: ubuntu-latest env: - LBPM_ZLIB_DIR: /home/runner/extlib/zlib - LBPM_HDF5_DIR: /home/runner/extlib/hdf5 - LBPM_SILO_DIR: /home/runner/extlib/silo - MPI_DIR: /home/runner/.openmpi - - steps: - - name: download dependencies - run: | - echo $LBPM_ZLIB_DIR - echo $LBPM_HDF5_DIR - echo $LBPM_SILO_DIR - echo $GITHUB_PATH - echo $GITHUB_WORKSPACE - - sudo apt-get update -y - - wget https://bitbucket.org/AdvancedMultiPhysics/tpl-builder/downloads/Silo-4.10.3RC.modified.tar.gz - wget https://www.zlib.net/fossils/zlib-1.2.11.tar.gz - wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.12/src/hdf5-1.8.12.tar.gz - #wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.10/src/hdf5-1.8.10.tar.gz - - tar -xzvf zlib-1.2.11.tar.gz - tar -xzvf hdf5-1.8.12.tar.gz - tar -xzvf Silo-4.10.3RC.modified.tar.gz - - - - name: check out commit - uses: actions/checkout@v2 - with: - path: LBPM - + MPI_VERSION: "4.1.8" + ZLIB_VERSION: "1.3.2" + HDF5_VERSION: "1.14.6" + MPI_DIR: ${{ github.workspace }}/deps/openmpi + LBPM_ZLIB_DIR: ${{ github.workspace }}/deps/zlib + LBPM_HDF5_DIR: ${{ github.workspace }}/hdf5 - - name: install-openmpi - run: | - wget https://download.open-mpi.org/release/open-mpi/v3.1/openmpi-3.1.2.tar.gz - tar -xvf ./openmpi-3.1.2.tar.gz - ./openmpi-3.1.2/configure --prefix="$HOME/.openmpi" - make -j - sudo make install - echo "$HOME/.openmpi/bin" >> $GITHUB_PATH - - + steps: + - uses: actions/checkout@v4 + with: + path: LBPM - - name: install zlib dependencies - run: | - cd zlib-1.2.11 - ./configure --prefix=$LBPM_ZLIB_DIR - make - sudo make install - cd .. + - name: Install build tools + run: | + sudo apt-get update -y + sudo apt-get install -y build-essential cmake wget tar + - name: Cache compiled deps (OpenMPI, zlib, HDF5) + id: cache-deps + uses: actions/cache@v4 + with: + path: | + deps/openmpi + deps/zlib + deps/hdf5 + key: ${{ runner.os }}-lbpm-deps-mpi${{ env.MPI_VERSION }}-zlib${{ env.ZLIB_VERSION }}-hdf5${{ env.HDF5_VERSION }} + restore-keys: | + ${{ runner.os }}-lbpm-deps- - - name: install hdf5 dependencies - run: | - cd hdf5-1.8.12 - CC=/home/runner/.openmpi/bin/mpicc CXX=/home/runner/.openmpi/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \ - ./configure --prefix=$LBPM_HDF5_DIR --enable-parallel --enable-shared --with-zlib=$LBPM_ZLIB_DIR - make - sudo make install - cd .. + - name: Build OpenMPI + if: steps.cache-deps.outputs.cache-hit != 'true' + run: | + wget -q https://download.open-mpi.org/release/open-mpi/v4.1/openmpi-${MPI_VERSION}.tar.gz + tar -xzf openmpi-${MPI_VERSION}.tar.gz + cd openmpi-${MPI_VERSION} + ./configure --prefix="${MPI_DIR}" + make -j + make install + - name: Build zlib + if: steps.cache-deps.outputs.cache-hit != 'true' + run: | + wget -q https://zlib.net/zlib-${ZLIB_VERSION}.tar.gz + tar -xzf zlib-${ZLIB_VERSION}.tar.gz + cd zlib-${ZLIB_VERSION} + ./configure --prefix="${LBPM_ZLIB_DIR}" + make -j + make install - - name: install silo dependencies - run: | - cd Silo-4.10.3RC.modified - CC=$MPI_DIR/bin/mpicc CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \ - ./configure --prefix=$LBPM_SILO_DIR -with-hdf5="$LBPM_HDF5_DIR/include,$LBPM_HDF5_DIR/lib" --enable-static - make - sudo make install - cd .. - - - - name: configure cmake - run: | - mkdir build - cd build - rm -rf CMake* - cmake \ - -D CMAKE_BUILD_TYPE:STRING=Release \ - -D CMAKE_C_COMPILER:PATH=$MPI_DIR/bin/mpicc \ - -D CMAKE_CXX_COMPILER:PATH=$MPI_DIR/bin/mpicxx \ - -D MPI_CXX_COMPILER=$MPI_DIR/bin/mpicxx \ - -D CMAKE_C_FLAGS="-fPIC" \ - -D CMAKE_CXX_FLAGS="-fPIC" \ - -D CMAKE_CXX_STD=14 \ - -D TEST_MAX_PROCS=1 \ - -D USE_TIMER=0 \ - -D TIMER_DIRECTORY=$LBPM_TIMER_DIR \ - -D USE_NETCDF=0 \ - -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR \ - -D USE_SILO=0 \ - -D HDF5_DIRECTORY=$LBPM_HDF5_DIR \ - -D SILO_DIRECTORY=$LBPM_SILO_DIR \ - -D USE_CUDA=0 \ - $GITHUB_WORKSPACE/LBPM + - name: Build HDF5 (parallel) + if: steps.cache-deps.outputs.cache-hit != 'true' + run: | + wget -q https://sourceforge.net/projects/hdf5.mirror/files/hdf5_${HDF5_VERSION}/hdf5-${HDF5_VERSION}.tar.gz/download -O hdf5-${HDF5_VERSION}.tar.gz + tar -xzf hdf5-${HDF5_VERSION}.tar.gz + cd hdf5-${HDF5_VERSION} + CC="${MPI_DIR}/bin/mpicc" CXX="${MPI_DIR}/bin/mpicxx" CXXFLAGS="-fPIC -O3 -std=c++14" \ + ./configure --prefix="${LBPM_HDF5_DIR}" --enable-parallel --enable-shared --with-zlib="${LBPM_ZLIB_DIR}" + make -j + make install + + # Ensure the cached/restored MPI is on PATH + - name: Add MPI to PATH + run: echo "${MPI_DIR}/bin" >> $GITHUB_PATH - + - name: Configure CMake + run: | + cmake -S LBPM -B build \ + -D USE_SILO=0 \ + -D USE_TIMER=0 \ + -D USE_CUDA=0 \ + -D CMAKE_BUILD_TYPE=Release \ + -D CMAKE_C_COMPILER:PATH="${MPI_DIR}/bin/mpicc" \ + -D CMAKE_CXX_COMPILER:PATH="${MPI_DIR}/bin/mpicxx" \ + -D MPI_CXX_COMPILER="${MPI_DIR}/bin/mpicxx" \ + -D CMAKE_C_FLAGS="-fPIC" \ + -D CMAKE_CXX_FLAGS="-fPIC" \ + -D CMAKE_CXX_STD=14 \ + -D TEST_MAX_PROCS=1 \ + -D HDF5_DIRECTORY="${LBPM_HDF5_DIR}" - - name: build and make - run: | - cd build - make - sudo make install - cd .. - - - name: tests - run: | - cd build - ctest + - name: build and make + run: | + cd build + make install + + - name: tests + run: | + ctest --test-dir build --output-on-failure diff --git a/.github/workflows/test_install_openmpi.yml b/.github/workflows/test_install_openmpi.yml deleted file mode 100644 index b193a6e6..00000000 --- a/.github/workflows/test_install_openmpi.yml +++ /dev/null @@ -1,41 +0,0 @@ -name: Install OpenMPI test - -on: - push: - branches: [ master ] - pull_request: - branches: [ master ] - -jobs: - install-openmpi: - - runs-on: ubuntu-latest - - steps: - - name: check path - run: | - echo $PATH - echo $GITHUB_PATH - cmake --version - - - name: download-openmpi - run: wget https://download.open-mpi.org/release/open-mpi/v4.0/openmpi-4.0.2.tar.gz - - - name: extract-openmpi - run: tar -xvf ./openmpi-4.0.2.tar.gz - - - name: configure-openmpi - run: ./openmpi-4.0.2/configure --prefix="/home/${USER}/.openmpi" - - - name: install-openmpi - run: | - make -j - sudo make install - - - name: setting path - run: | - echo "/home/${USER}/.openmpi/bin" >> $GITHUB_PATH - #echo "/home/${USER}/.openmpi/bin" >> $PATH - - name: checking version - run: mpirun --version - diff --git a/IO/Mesh.cpp b/IO/Mesh.cpp index 5e5eb96a..0eb9283f 100644 --- a/IO/Mesh.cpp +++ b/IO/Mesh.cpp @@ -598,6 +598,8 @@ std::string getString( FileFormat type ) return "silo"; else if ( type == FileFormat::HDF5 ) return "hdf5"; + else if ( type == FileFormat::VTK ) + return "vtk"; else ERROR( "Invalid type" ); return ""; diff --git a/IO/Mesh.h b/IO/Mesh.h index 9e5f32e6..2ada354c 100644 --- a/IO/Mesh.h +++ b/IO/Mesh.h @@ -24,7 +24,7 @@ enum class VariableType { }; enum class DataType { Double, Float, Int, Null }; enum class MeshType { PointMesh, SurfaceMesh, VolumeMesh, Unknown }; -enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 }; +enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 , VTK}; //! Convert enums to/from strings (more future-proof than static_cast) diff --git a/IO/Reader.cpp b/IO/Reader.cpp index 014e1e7b..410b78a4 100644 --- a/IO/Reader.cpp +++ b/IO/Reader.cpp @@ -5,6 +5,7 @@ #include "IO/MeshDatabase.h" #include "IO/silo.h" #include "common/Utilities.h" +#include "IO/silo.h" #include #include diff --git a/IO/VtiWriter.cpp b/IO/VtiWriter.cpp new file mode 100755 index 00000000..1e6223e0 --- /dev/null +++ b/IO/VtiWriter.cpp @@ -0,0 +1,117 @@ +#include "IO/HDF5_IO.h" +#include "IO/IOHelpers.h" +#include "IO/MeshDatabase.h" +#include "IO/Writer.h" +#include "IO/silo.h" +#include "IO/xmlvtk.h" +#include "common/MPI.h" +#include "common/Utilities.h" + +#include +#include +#include +#include +#include + +static void writeVti( + const std::string &fullpath, const IO::MeshDataStruct &meshData) +{ + const IO::DomainMesh &mesh = dynamic_cast( *meshData.mesh ); + RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + + VTIWriter vti = VTIWriter(std::string(fullpath)); + vti.setWholeExtent( info.ix * mesh.nx, info.jy * mesh.ny , info.kz * mesh.nz, + ( info.ix + 1 ) * mesh.nx, ( info.jy + 1 ) * mesh.ny, ( info.kz + 1 ) * mesh.nz); + + vti.setSpacing( 1.0 , 1.0 , 1.0 ); + vti.setOrigin(0,0,0); + vti.setCompress(); + + for ( size_t i = 0; i < meshData.vars.size(); i++ ) + { + const auto &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + vti.addCellData( var.name , "Float64" , "binary" , var.dim , (unsigned char*) var.data.begin() ); + } else if ( var.precision == IO::DataType::Float ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + vti.addCellData( var.name , "Float32" , "binary" , var.dim , (unsigned char*) data2.begin()); + } else if ( var.precision == IO::DataType::Int ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + vti.addCellData( var.name , "Int32" , "binary" , var.dim , (unsigned char*) var.data.begin() ); + } else { + ERROR( "Unsupported format" ); + } + } + + vti.write(); +} + +void writeVtiSummary( + const std::vector &meshes_written,const IO::MeshDataStruct &meshData, const std::string &filename ) +{ + const IO::DomainMesh &mesh = dynamic_cast( *meshData.mesh ); + RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + PVTIWriter pvti = PVTIWriter( filename ); + int rank = 0; + for ( const auto &data : meshes_written ) + { + for ( const auto &tmp : data.domains ) + { + RankInfoStruct info( rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + char filename[100]; + sprintf( filename, "%05i.vti", rank ); + + VTIWriter vti = VTIWriter( filename ); + vti.setWholeExtent( info.ix * mesh.nx, info.jy * mesh.ny , info.kz * mesh.nz, + ( info.ix + 1 ) * mesh.nx, ( info.jy + 1 ) * mesh.ny, ( info.kz + 1 ) * mesh.nz); + vti.setSpacing( 1.0 , 1.0 , 1.0 ); + vti.setOrigin(0,0,0); + vti.setCompress(); + + for ( size_t i = 0; i < meshData.vars.size(); i++ ) + { + const auto &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + vti.addCellData( var.name , "Float64" , "binary" , var.dim , nullptr ); + } else if ( var.precision == IO::DataType::Float ) { + vti.addCellData( var.name , "Float32" , "binary" , var.dim , nullptr ); + } else if ( var.precision == IO::DataType::Int ) { + vti.addCellData( var.name , "Int32" , "binary" , var.dim , nullptr ); + } else { + ERROR( "Unsupported format" ); + } + } + + pvti.addVTIWriter(vti); + rank++; + + } + } + pvti.write(); +} + +std::vector writeMeshesVti( const std::vector &meshData, + const std::string &path, int rank ) +{ + std::vector meshes_written; + char filename[100], fullpath[200]; + sprintf( filename, "%05i.vti", rank ); + sprintf( fullpath, "%s/%s", path.c_str(), filename ); + + for ( size_t i = 0; i < meshData.size(); i++ ) { +// auto mesh = meshData[i].mesh; + auto database = getDatabase( fullpath , meshData[i], IO::FileFormat::VTK, rank ); + + if ( database.meshClass == "DomainMesh" ) { + writeVti( fullpath, meshData[i] ); + } else { + ERROR( "Unknown mesh class or not implemented for vtk/vti output" ); + } + + meshes_written.push_back( database ); + } + return meshes_written; +} + diff --git a/IO/Writer.cpp b/IO/Writer.cpp index 0fb8d135..1ed3cdfa 100644 --- a/IO/Writer.cpp +++ b/IO/Writer.cpp @@ -15,7 +15,7 @@ #include -enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN }; +enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN, VTK }; /**************************************************** @@ -26,7 +26,10 @@ std::vector writeMeshesSilo( void writeSiloSummary( const std::vector &, const std::string & ); std::vector writeMeshesHDF5( const std::vector &, const std::string &, IO::FileFormat, int, Xdmf & ); - +std::vector writeMeshesVti( const std::vector &meshData, + const std::string &path, int rank ); +void writeVtiSummary( + const std::vector &meshes_written,const IO::MeshDataStruct &meshData, const std::string &filename ); /**************************************************** * Recursively create the subdirectory * @@ -90,6 +93,8 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap global_IO_format = Format::SILO; else if ( format == "hdf5" ) global_IO_format = Format::HDF5; + else if ( format == "vtk" ) + global_IO_format = Format::VTK; else ERROR( "Unknown format" ); int rank = Utilities::MPI( MPI_COMM_WORLD ).getRank(); @@ -100,9 +105,13 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap filename = global_IO_path + "/summary.LBM"; else if ( global_IO_format == Format::SILO || global_IO_format == Format::HDF5 ) filename = global_IO_path + "/LBM.visit"; + else if ( global_IO_format == Format::VTK) + filename = global_IO_path + "/LBM.pvd"; else ERROR( "Unknown format" ); auto fid = fopen( filename.c_str(), "wb" ); + if ( global_IO_format == Format::VTK) + fprintf( fid, "\n\n \n \n\n" ); fclose( fid ); } } @@ -278,13 +287,12 @@ static std::vector writeMeshesNewFormat( return meshes_written; } - /**************************************************** * Write the mesh data * ****************************************************/ void IO::writeData( const std::string &subdir, const std::vector &meshData, - const Utilities::MPI &comm ) -{ + const Utilities::MPI &comm , int timestep ) +{ if ( global_IO_path.empty() ) IO::initialize(); PROFILE_START( "writeData" ); @@ -310,6 +318,9 @@ void IO::writeData( const std::string &subdir, const std::vector\n \n\n", timestep, subdir.c_str()); + + char *pos = strstr(buffer, " \n\n"); + pos[0] = '\0'; + + fid = fopen( filename.c_str(), "w" ); + fputs(buffer,fid); + fputs(newrow,fid); + fclose(fid); + + free(buffer); + } else { ERROR( "Unknown format" ); } diff --git a/IO/Writer.h b/IO/Writer.h index 3844f3b2..e1d2f8f4 100644 --- a/IO/Writer.h +++ b/IO/Writer.h @@ -38,7 +38,7 @@ void initialize( * @param[in] comm The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof) */ void writeData( const std::string &subdir, const std::vector &meshData, - const Utilities::MPI &comm ); + const Utilities::MPI &comm , int timestep = -1); /*! @@ -48,12 +48,14 @@ void writeData( const std::string &subdir, const std::vector * @param[in] meshData The data to write * @param[in] comm The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof) */ + + inline void writeData( int timestep, const std::vector &meshData, const Utilities::MPI &comm ) { char subdir[100]; sprintf( subdir, "vis%03i", timestep ); - writeData( subdir, meshData, comm ); + writeData( subdir, meshData, comm , timestep); } diff --git a/IO/xmlvtk.cpp b/IO/xmlvtk.cpp new file mode 100644 index 00000000..aceb4e30 --- /dev/null +++ b/IO/xmlvtk.cpp @@ -0,0 +1,559 @@ +/* + * Copyright (c) 2025 Diogo Nardelli Siebert + * + * Licensed under either of + * - Apache License, Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0) + * - GNU General Public License, Version 3.0 or later (https://www.gnu.org/licenses/gpl-3.0.html) + * + * SPDX-License-Identifier: (Apache-2.0 OR GPL-3.0-or-later) + */ + +#include "xmlvtk.h" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +/** @brief Indicates whether data compression is used by default **/ +bool Element::compress = false; +/** @brief Default cache size used during data processing. */ +headerType cacheSize = 10000; + +/** + * @brief Formats a string using printf-style syntax. + * @tparam Args Variadic arguments for formatting. + * @param fmt Format string. + * @return Formatted string. + */ +template +std::string format_string(const char* fmt, Args&&... args) +{ + int size = std::snprintf(nullptr, 0, fmt, std::forward(args)...); + if (size < 0) { + throw std::runtime_error("format_string: snprintf error"); + } + std::vector buf(size + 1); + int size2 = std::snprintf(buf.data(), buf.size(), fmt, std::forward(args)...); + if (size2 < 0) { + throw std::runtime_error("format_string: snprintf error"); + } + + return std::string(buf.data(), buf.data() + size2); +} + +/** + * @brief Encodes binary input to Base64 string representation. + * @param input Input binary data. + * @param len Length of input data. + * @return Encoded Base64 string. + */ +std::string spc_base64_encode(const unsigned char* input, size_t len) +{ + static const char table[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/"; + std::string output; + output.reserve((len / 3 + (len % 3 != 0)) * 4); + + for (size_t i = 0; i < len; i += 3) { + uint32_t n = (input[i] << 16) | (i + 1 < len ? input[i + 1] << 8 : 0) | (i + 2 < len ? input[i + 2] : 0); // AAAAAABB BBBBCCCC CCDDDDDD + output.push_back(table[(n >> 18) & 0x3F]); + output.push_back(table[(n >> 12) & 0x3F]); + output.push_back(i + 1 < len ? table[(n >> 6) & 0x3F] : '='); + output.push_back(i + 2 < len ? table[n & 0x3F] : '='); + } + return output; +} + +std::string CellData::header() +{ + return "\n"; +} + +std::string CellData::footer() +{ + return "\n"; +} + +std::string PointData::header() +{ + return "\n"; +} + +std::string PointData::footer() +{ + return "\n"; +} + +std::string AppendData::header() +{ + return "\n"; +} + +std::string AppendData::footer() +{ + return "\n"; +} + +AppendData::AppendData() +{ + this -> totalSize = 0; +} + +unsigned int AppendData::addData(unsigned char* pointer, unsigned int size) +{ + int offset = this->totalSize; + this->sizeList.push_back(size); + this->pointerList.push_back(pointer); + this->totalSize += (size+4); + return offset; +} + +DataArray::DataArray(const std::string& name_,const std::string& type_, const std::string& format_, int components_, unsigned char* pointer_, uint64_t points_) + : components(components_), name(name_), type(type_), format(format_), pointer(pointer_), points(points_) + +{ + if ( (this -> type == "Int8") || (this -> type == "UInt8") ) this -> typeSize = 1; + else if ( (this -> type == "Int16") || (this -> type == "UInt16") ) this -> typeSize = 2; + else if ( (this -> type == "Int32") || (this -> type == "UInt32") || (this -> type == "Float32") ) this -> typeSize = 4; + else if ( (this -> type == "Int64") || (this -> type == "UInt64") || (this -> type == "Float64") ) this -> typeSize = 8; + this -> dataSize = components * points; +} + +std::string DataArray::header() +{ + std::ostringstream stringStream; + stringStream << format_string("", offset).c_str() : ">"); + return stringStream.str(); +} + +std::string DataArray::footer() +{ + return "\n\n"; +} + +void VTIWriter::setCompress() +{ + this->compress = true; +} + +PVTIWriter::PVTIWriter(const std::string& filename_) + : filename(filename_), + originX(0), originY(0), originZ(0), + wholeMinX(0), wholeMinY(0), wholeMinZ(0), wholeMaxX(0), wholeMaxY(0), wholeMaxZ(0), pieceCounter(0) +{ +} + +VTIWriter::VTIWriter(const std::string& filename_) + : filename(filename_) , vtkVersion("0.1"), + originX(0), originY(0), originZ(0) +{ + this -> compress = false; + if (sizeof(headerType) == 8) this -> headerTypeName = "UInt64"; + else if (sizeof(headerType) == 4) this -> headerTypeName = "UInt32"; + this -> byteOrder = "LittleEndian"; +} + +void VTIWriter::addPointData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer) +{ + std::unique_ptr data = std::make_unique( name , type , format , components , (unsigned char*) pointer , pd.points ); + data -> compress = this -> compress; + + if (format == "appended") + { + data -> offset = ad.addData( data->pointer , data->dataSize * data->typeSize ); + ad.compress = this->compress; + } + + pd.addChild( std::move(data) ); +} + +void VTIWriter::addCellData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer) +{ + std::unique_ptr data = std::make_unique( name , type , format , components , (unsigned char*) pointer , cd.cells ); + data -> compress = this -> compress; + + if (format == "appended") + { + data -> offset = ad.addData( data->pointer , data->dataSize * data->typeSize ); + ad.compress = this->compress; + } + + cd.addChild( std::move(data) ); +} + +std::string VTIWriter::header() +{ + std::ostringstream stringStream; + stringStream << "" << endl; + stringStream << format_string("", "1.0", byteOrder.c_str() , headerTypeName.c_str() , compress ? "compressor=\"vtkZLibDataCompressor\"" : "" ) << endl; + stringStream << format_string("", spaceX, spaceY, spaceZ) << endl; + stringStream << format_string("" , pieceMinX , pieceMaxX , pieceMinY , pieceMaxY, pieceMinZ , pieceMaxZ) << endl; + return stringStream.str(); +} + +std::string PVTIWriter::header() +{ + std::ostringstream stringStream; + stringStream << "" << endl; + stringStream << "" << endl; + stringStream << format_string("", spaceX, spaceY, spaceZ) << endl; + + return stringStream.str(); +} + +std::string VTIWriter::footer() +{ + return "\n"; +} + +std::string PVTIWriter::footer() +{ + return "\n\n"; +} + +void VTIWriter::setWholeExtent(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_) +{ + wholeMinX = minX_; + wholeMinY = minY_; + wholeMinZ = minZ_; + wholeMaxX = maxX_; + wholeMaxY = maxY_; + wholeMaxZ = maxZ_; + setPiece(minX_,minY_,minZ_,maxX_,maxY_, maxZ_); +} + +std::ostream& operator<<(std::ostream& os, AppendData& obj) +{ + if (obj.totalSize > 0) + { + os << obj.header(); + os << "_"; + for (int k = 0; k < obj.pointerList.size(); k++) + { + if (obj.compress == false) + { + os.write( (char*) &obj.sizeList[k], sizeof(headerType) ); + os.write( (char*) obj.pointerList[k], obj.sizeList[k] ); + } + else + { + streampos beginPos = os.tellp(); + + headerType totalByteSize = obj.sizeList[k]; + headerType numberOfBlocks = totalByteSize / cacheSize + (totalByteSize % cacheSize != 0 ); + + size_t infoSize = sizeof(headerType) * (3 + numberOfBlocks); + std::vectorcompressedInfo(infoSize ); + + compressedInfo[0] = numberOfBlocks; + compressedInfo[1] = (numberOfBlocks > 1) ? cacheSize : totalByteSize; + compressedInfo[2] = (totalByteSize % cacheSize == 0) ? cacheSize : totalByteSize % cacheSize ; + + os.write((char*) compressedInfo.data() , infoSize ); + + for (int n = 1; n <= numberOfBlocks; n++) + { + int numberOfBytesInBlock = (n < numberOfBlocks) ? compressedInfo[1] : compressedInfo[2]; + + uLongf compressedLength = compressBound(numberOfBytesInBlock); + std::vectorcompressedData (compressedLength); + + const Bytef* src = reinterpret_cast( obj.pointerList[k] + (n - 1) * cacheSize ); + int ret = compress( (Bytef *) (compressedData.data() ), &compressedLength, src, numberOfBytesInBlock ); + if (ret != Z_OK) + { + std::cerr << "Error compressing data , code =" << ret << std::endl; + throw std::runtime_error("Failed compressing data using zlib in AppendData"); + } + compressedInfo[n+2] = compressedLength; + + os.write( reinterpret_cast( compressedData.data() ) , compressedLength ); + } + + streampos endPos = os.tellp(); + + os.seekp( beginPos ); + os.write( (char*) compressedInfo.data() , infoSize ); + os.seekp( endPos ); + } + } + os << endl; + os << obj.footer(); + + } + return os; +} + +std::ostream& operator<<(std::ostream& os, DataArray& obj) +{ + os << obj.header() << endl; + + headerType totalByteSize = static_cast( (obj.dataSize) * (obj.typeSize) ); + headerType count = min( 12 - sizeof(headerType) , totalByteSize ); + + if (obj.format == "binary") + { + if (obj.compress == false) + { + os << spc_base64_encode( (unsigned char *) &totalByteSize, sizeof(headerType) ); + count = 0; + + int leftOverSize = 0; + unsigned char leftOverBuffer[3]; + + for (; count < totalByteSize ; count += cacheSize) + { + unsigned char* pointer = obj.pointer + count; + size_t size = min( cacheSize, totalByteSize - count); + if (leftOverSize > 0) + { + for (; leftOverSize < 3; leftOverSize++) + { + if (size > 0) + { + leftOverBuffer[leftOverSize] = *(pointer++); + size--; + } + else break; + } + + os << spc_base64_encode( leftOverBuffer , leftOverSize ); + leftOverSize = 0; + } + + leftOverSize = (size % 3); + size -= leftOverSize; + for (int i = 0; i < leftOverSize; i++) + { + leftOverBuffer[i] = pointer[size + i]; + } + + os << spc_base64_encode( pointer , size ); + + } + + if (leftOverSize > 0) + { + os << spc_base64_encode( leftOverBuffer , leftOverSize ); + leftOverSize = 0; + } + } + else + { + streampos beginPos = os.tellp(); + + std::vectorcompressedData (cacheSize + 3); + int leftOver = 0; + + headerType numberOfBlocks = totalByteSize / cacheSize + (totalByteSize % cacheSize != 0 ); + size_t infoSize = sizeof(headerType) * (3 + numberOfBlocks); + std::vectorcompressedInfo(infoSize ); + + compressedInfo[0] = numberOfBlocks; + compressedInfo[1] = (numberOfBlocks > 1) ? cacheSize : totalByteSize; + compressedInfo[2] = (totalByteSize % cacheSize == 0) ? cacheSize : totalByteSize % cacheSize ; + + os << spc_base64_encode( reinterpret_cast(compressedInfo.data()) , infoSize ); + + for (int n = 1; n <= numberOfBlocks; n++) + { + uLongf compressedLength = cacheSize; + int numberOfBytesInBlock = (n < numberOfBlocks) ? compressedInfo[1] : compressedInfo[2]; + int ret = compress( reinterpret_cast(compressedData.data() + 3) , &compressedLength, obj.pointer + (n-1) * cacheSize, numberOfBytesInBlock ); + if (ret != Z_OK) + { + std::cerr << "Error compressing data , code =" << ret << std::endl; + throw std::runtime_error("Failed compressing data using zlib in AppendData"); + } + compressedInfo[n+2] = compressedLength; + + int encodeSize = 3* ( (compressedLength + leftOver)/3 ); + os << spc_base64_encode( compressedData.data() + 3 - leftOver , encodeSize ); + + int newLeftOver = compressedLength + leftOver - encodeSize ; + if (newLeftOver > 0) + { + memcpy( compressedData.data() + 3 - newLeftOver, compressedData.data() + 3 - leftOver + encodeSize, newLeftOver); + } + leftOver = newLeftOver; + } + + if (leftOver > 0) + { + os << spc_base64_encode(compressedData.data() + 3 - leftOver , leftOver ); + } + + streampos endPos = os.tellp(); + os.seekp( beginPos ); + + os << spc_base64_encode( reinterpret_cast(compressedInfo.data()) , infoSize ); + os.seekp( endPos ); + } + + os << obj.footer(); + } + + return os; +} + +void VTIWriter::setPiece(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_) +{ + pieceMinX = minX_; + pieceMinY = minY_; + pieceMinZ = minZ_; + pieceMaxX = maxX_; + pieceMaxY = maxY_; + pieceMaxZ = maxZ_; + sizeX = maxX_ - minX_ + 1; + sizeY = maxY_ - minY_ + 1; + sizeZ = maxZ_ - minZ_ + 1; + pd.points = sizeX * sizeY * sizeZ; + cd.cells = (sizeX-1) * (sizeY-1) * (sizeZ-1); +} + +void VTIWriter::setOrigin(double x_,double y_, double z_) +{ + originX = x_; + originY = y_; + originZ = z_; +} + +void VTIWriter::setSpacing(double sx_,double sy_, double sz_) +{ + spaceX = sx_; + spaceY = sy_; + spaceZ = sz_; +} + +void PVTIWriter::write() +{ + std::ofstream file; + file.open(this -> filename); + file << header(); + + for (int i = 0; i< cellDataName.size(); i++) + { + if (i==0) file << "" << endl; + file << format_string( "" ,cellDataName[i].c_str() ,cellDataComponents[i] ,cellDataType[i].c_str() ) << endl; + if (i== cellDataName.size()-1) file << "" << endl; + } + + for (int i = 0; i< pointDataName.size(); i++) + { + if (i==0) file << "" << endl; + file << format_string( "" ,pointDataName[i].c_str() ,pointDataComponents[i] ,pointDataType[i].c_str() ) << endl; + if (i== pointDataName.size()-1) file << "" << endl; + } + + + for (int i = 0; i< pieceFilename.size(); i++) + { + file << format_string("", pieceMinX[i], pieceMaxX[i], pieceMinY[i] , pieceMaxY[i] , pieceMinZ[i] , pieceMaxZ[i] , pieceFilename[i].c_str() ) << endl; + } + + file << footer(); + file.close(); +} + +void VTIWriter::write() +{ + file.open(this -> filename); + file << header(); + + if (pd.sizeChild() > 0) + { + file << pd.header(); + for (int n = 0; n < pd.sizeChild() ; n++) + { + DataArray* array = (DataArray*) pd.getChild(n); + file << *array; + } + file << pd.footer(); + } + + if (cd.sizeChild() > 0) + { + file << cd.header(); + for (int n = 0; n < cd.sizeChild() ; n++) + { + DataArray* array = (DataArray*) cd.getChild(n); + file << *array; + } + file << cd.footer(); + } + + file << "" << endl; + file << "" << endl; + + file << ad; + file << footer(); + file.close(); +} + +void PVTIWriter::addVTIWriter(VTIWriter& write) +{ + pieceFilename.push_back( write.filename ); + + pieceMaxX.push_back( write.pieceMaxX ); + pieceMaxY.push_back( write.pieceMaxY ); + pieceMaxZ.push_back( write.pieceMaxZ ); + pieceMinX.push_back( write.pieceMinX ); + pieceMinY.push_back( write.pieceMinY ); + pieceMinZ.push_back( write.pieceMinZ ); + + if (pieceFilename.size() == 1) + { + for (int n = 0; n < write.pd.sizeChild() ; n++) + { + DataArray* data = (DataArray*) write.pd.getChild(n); + pointDataName.push_back( data -> name ); + pointDataType.push_back( data -> type); + pointDataComponents.push_back( data -> components ); + } + + for (int n = 0; n < write.cd.sizeChild() ; n++) + { + DataArray* data = (DataArray*) write.cd.getChild(n); + cellDataName.push_back( data -> name ); + cellDataType.push_back( data -> type); + cellDataComponents.push_back( data -> components ); + } + + originX = write.originX; + originY = write.originX; + originZ = write.originX; + + spaceX = write.spaceX; + spaceY = write.spaceY; + spaceZ = write.spaceZ; + + wholeMaxX = write.wholeMaxX; + wholeMaxY = write.wholeMaxY; + wholeMaxZ = write.wholeMaxZ; + wholeMinX = write.wholeMinX; + wholeMinY = write.wholeMinY; + wholeMinZ = write.wholeMinZ; + } + else + { + wholeMaxX = max(wholeMaxX,write.wholeMaxX); + wholeMaxY = max(wholeMaxY,write.wholeMaxY); + wholeMaxZ = max(wholeMaxZ,write.wholeMaxZ); + wholeMinX = min(wholeMinX,write.wholeMinX); + wholeMinY = min(wholeMinY,write.wholeMinY); + wholeMinZ = min(wholeMinZ,write.wholeMinZ); + } + + +} diff --git a/IO/xmlvtk.h b/IO/xmlvtk.h new file mode 100644 index 00000000..9538eda4 --- /dev/null +++ b/IO/xmlvtk.h @@ -0,0 +1,185 @@ +/* + * Copyright (c) 2025 Diogo Nardelli Siebert + * + * Licensed under either of + * - Apache License, Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0) + * - GNU General Public License, Version 3.0 or later (https://www.gnu.org/licenses/gpl-3.0.html) + * + * SPDX-License-Identifier: (Apache-2.0 OR GPL-3.0-or-later) + */ + +#ifndef __XMLVTK_H_INCLUDED__ +#define __XMLVTK_H_INCLUDED__ + +#include +#include +#include +#include +#include + +#define headerType u_int64_t + + +/** + * @brief Encodes binary data into Base64 format. + * @param input Pointer to the input data. + * @param len Length of the input data. + * @return Base64-encoded string. + */ +std::string spc_base64_encode(const unsigned char* input, size_t len); + +/** + * @class Element + * @brief Represents a hierarchical XML-like element supporting compression and VTK serialization. + */ +class Element +{ + public: + static bool compress; +/** @brief Returns the parent element. */ + Element* getParent() { return parent ;} +/** @brief Returns the nth child element. *//** ... */ + Element* getChild(int n) { return this->child[n].get(); } +/** @brief Adds a new child element and sets its parent. */ + void addChild(std::unique_ptr e) { + e -> parent = this; + child.push_back(std::move( e ) ); + } + +/** @brief Returns the number of child elements. */ + int sizeChild( ) { return child.size(); } +/** @brief Generates a header string for VTK output. */ + virtual std::string header() = 0; +/** @brief Generates a footer string for VTK output. */ + virtual std::string footer() = 0; + private: + Element* parent; + std::vector< std::unique_ptr > child; +}; + +class PointData: public Element +{ + public: + std::string header(); + std::string footer(); + unsigned int points; +}; + +class CellData: public Element +{ + public: + std::string header(); + std::string footer(); + unsigned int cells; +}; + +class DataArray : public Element +{ + public: + DataArray(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer, uint64_t size); + unsigned char* pointer; + std::string header(); + std::string footer(); + void write( std::ofstream& file ); + friend std::ostream& operator<<(std::ostream& os, DataArray& obj); + uint64_t typeSize; + uint64_t dataSize ; + uint64_t offset; + uint64_t points; + std::string name; + std::string format; + std::string type; + int components; + private: + int mode = 0; + + +}; + +class AppendData : public Element +{ + public: + AppendData(); + std::string header(); + std::string footer(); + unsigned int addData(unsigned char* pointer, unsigned int size); + friend std::ostream& operator<<(std::ostream& os, AppendData& obj); + private: + std::vector pointerList; + std::vector sizeList; + headerType totalSize; +}; + +class VTIWriter : public Element +{ + public: + VTIWriter(const std::string& filename); + std::string footer(); + std::string header(); + + void write(); + void setWholeExtent(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_); + void setPiece(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_); + void setOrigin(double x_,double y_, double z_); + void setSpacing(double sx_,double sy_, double sz_); + void addPointData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer); + void addCellData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer); + void setCompress(); + + std::string filename; + std::string vtkVersion; /*!< The Vtk File Format Version of the file */ + std::string fileTitle; /*!< The title of the file (do not confuse with the name of the file) */ + std::string dataSetType; /*!< The type of geometry (grid) that data is associeted to (STRUCTURED GRID for LBM applications) */ + std::string headerTypeName; + + int64_t wholeMinX, wholeMinY, wholeMinZ; + int64_t wholeMaxX, wholeMaxY, wholeMaxZ; + int64_t pieceMinX, pieceMinY, pieceMinZ; + int64_t pieceMaxX, pieceMaxY, pieceMaxZ; + int64_t sizeX, sizeY, sizeZ; /*!< Lenght in pixels of the image in each axis */ + + /*!< Lenght in pixels of the image in each axis */ + double spaceX, spaceY, spaceZ; /*!< Ratio of the different axis */ + double originX,originY,originZ; /*!< Position of the origin of the image */ + + std::string byteOrder; /*!< Position of the origin of the image */ + std::ofstream file; + CellData cd; + PointData pd; + AppendData ad; +}; + +class PVTIWriter : public Element +{ + public: + PVTIWriter(const std::string& filename); + std::string filename; + std::string footer(); + std::string header(); + + double spaceX, spaceY, spaceZ; + double originX,originY,originZ; + + int pieceCounter; + + int64_t wholeMinX, wholeMinY, wholeMinZ; + int64_t wholeMaxX, wholeMaxY, wholeMaxZ; + + std::vector pieceMinX, pieceMinY, pieceMinZ; + std::vector pieceMaxX, pieceMaxY, pieceMaxZ; + std::vector pieceFilename; + + std::vector cellDataComponents; + std::vector cellDataName; + std::vector cellDataType; + + std::vector pointDataComponents; + std::vector pointDataName; + std::vector pointDataType; + + void addVTIWriter(VTIWriter& write); + void write(); +}; + + +#endif diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 728bd6ed..6a01812c 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -508,7 +508,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, // total Global is the number of nodes in the pore-space totalGlobal = Dm->Comm.sumReduce(count); - maxdistGlobal = Dm->Comm.sumReduce(maxdist); + maxdistGlobal = Dm->Comm.maxReduce(maxdist); double volume = double(nprocx * nprocy * nprocz) * double(nx - 2) * double(ny - 2) * double(nz - 2); double volume_fraction = totalGlobal / volume; @@ -683,7 +683,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, // total Global is the number of nodes in the pore-space totalGlobal = Dm->Comm.sumReduce(count); - maxdistGlobal = Dm->Comm.sumReduce(maxdist); + maxdistGlobal = Dm->Comm.maxReduce(maxdist); double volume = double(nprocx * nprocy * nprocz) * double(nx - 2) * double(ny - 2) * double(nz - 2); double volume_fraction = totalGlobal / volume; diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 79ff0834..e62b33ec 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -248,46 +248,53 @@ class IOWorkItem : public ThreadPool::WorkItemRet { fillHalo fillData(comm.comm, rank_info, n, {1, 1, 1}, 0, 1); + int varCounter = 0; if (vis_db->getWithDefault("save_phase_field", true)) { - ASSERT(visData[0].vars[0]->name == "phase"); - Array &PhaseData = visData[0].vars[0]->data; + ASSERT(visData[0].vars[varCounter]->name == "phase"); + Array &PhaseData = visData[0].vars[varCounter]->data; fillData.copy(Averages.Phi, PhaseData); + varCounter+=1; } if (vis_db->getWithDefault("save_pressure", false)) { - ASSERT(visData[0].vars[1]->name == "Pressure"); - Array &PressData = visData[0].vars[1]->data; + ASSERT(visData[0].vars[varCounter]->name == "Pressure"); + Array &PressData = visData[0].vars[varCounter]->data; fillData.copy(Averages.Pressure, PressData); + varCounter+=1; } if (vis_db->getWithDefault("save_velocity", false)) { - ASSERT(visData[0].vars[2]->name == "Velocity_x"); - ASSERT(visData[0].vars[3]->name == "Velocity_y"); - ASSERT(visData[0].vars[4]->name == "Velocity_z"); - Array &VelxData = visData[0].vars[2]->data; - Array &VelyData = visData[0].vars[3]->data; - Array &VelzData = visData[0].vars[4]->data; + ASSERT(visData[0].vars[varCounter]->name == "Velocity_x"); + ASSERT(visData[0].vars[varCounter+1]->name == "Velocity_y"); + ASSERT(visData[0].vars[varCounter+2]->name == "Velocity_z"); + Array &VelxData = visData[0].vars[varCounter]->data; + Array &VelyData = visData[0].vars[varCounter+1]->data; + Array &VelzData = visData[0].vars[varCounter+2]->data; fillData.copy(Averages.Vel_x, VelxData); fillData.copy(Averages.Vel_y, VelyData); fillData.copy(Averages.Vel_z, VelzData); + varCounter+=3; } if (vis_db->getWithDefault("save_dissipation", false)) { - ASSERT(visData[0].vars[5]->name == "ViscousDissipation"); - Array &ViscousDissipation = visData[0].vars[5]->data; + ASSERT(visData[0].vars[varCounter]->name == "ViscousDissipation"); + Array &ViscousDissipation = visData[0].vars[varCounter]->data; fillData.copy(Averages.Dissipation, ViscousDissipation); + varCounter+1; } if (vis_db->getWithDefault("save_distance", false)) { - ASSERT(visData[0].vars[6]->name == "SignDist"); - Array &SignData = visData[0].vars[6]->data; + ASSERT(visData[0].vars[varCounter]->name == "SignDist"); + Array &SignData = visData[0].vars[varCounter]->data; fillData.copy(Averages.SDs, SignData); + varCounter+=1; } if (vis_db->getWithDefault("save_connected_components", false)) { - ASSERT(visData[0].vars[7]->name == "BlobID"); - Array &BlobData = visData[0].vars[7]->data; + ASSERT(visData[0].vars[varCounter]->name == "BlobID"); + Array &BlobData = visData[0].vars[varCounter]->data; fillData.copy(Averages.morph_n->label, BlobData); + varCounter+=1; } if (vis_db->getWithDefault("write_silo", true)) @@ -591,7 +598,7 @@ runAnalysis::runAnalysis(std::shared_ptr input_db, //std::string format = "silo"; format = vis_db->getWithDefault("format", "silo"); - IO::initialize("", format, "false"); + IO::initialize("", format, false); // Create the MeshDataStruct d_meshData.resize(1); @@ -736,7 +743,7 @@ runAnalysis::runAnalysis(ScaLBL_ColorModel &ColorModel) format = vis_db->getWithDefault("format", "silo"); - IO::initialize("", format, "false"); + IO::initialize("", format, false); // Create the MeshDataStruct d_meshData.resize(1); diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 7bdc360a..521bda87 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -188,6 +188,8 @@ extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np); */ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np); +extern "C" void ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int Np, double Fx, double Fy, double Fz); + /** * \brief compute pressure from D3Q19 distribution * @param dist - D3Q19 distributions @@ -626,7 +628,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_MRT(int *neighborList, double *dist, * @param Np - size of local sub-domain (derived from Domain structure) */ extern "C" void ScaLBL_D3Q19_AAeven_Color( - int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); @@ -658,7 +660,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( */ extern "C" void ScaLBL_D3Q19_AAodd_Color( int *NeighborList, int *Map, double *dist, double *Aq, double *Bq, - double *Den, double *Phi, double *Vel, double rhoA, double rhoB, + double *Den, double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 526bcafa..15d84f78 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -17,6 +17,7 @@ #include #define STOKES +#define f64_eps 1E-12 extern "C" void ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, @@ -1429,7 +1430,7 @@ extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, // double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, // double Fx, double Fy, double Fz, int start, int finish, int Np){ extern "C" void ScaLBL_D3Q19_AAeven_Color( - int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int* NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np) { @@ -1446,6 +1447,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( double C, nx, ny, nz; //color gradient magnitude and direction double ux, uy, uz; double phi, tau, rho0, rlx_setA, rlx_setB; + double nspx, nspy, nspz; // Rock Fluid interface normal vector const double mrt_V1 = 0.05263157894736842; const double mrt_V2 = 0.012531328320802; @@ -1472,83 +1474,243 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = tauA + 0.5 * (1.0 - phi) * (tauB - tauA); + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); - // Get the 1D index based on regular data layout - ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ - nn = ijk - 1; // neighbor index (get convention) - m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ - nn = ijk + 1; // neighbor index (get convention) - m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ - nn = ijk - strideY; // neighbor index (get convention) - m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ - nn = ijk + strideY; // neighbor index (get convention) - m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ - nn = ijk - strideZ; // neighbor index (get convention) - m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ - nn = ijk + strideZ; // neighbor index (get convention) - m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ - nn = ijk - strideY - 1; // neighbor index (get convention) - m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ - nn = ijk + strideY + 1; // neighbor index (get convention) - m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ - nn = ijk + strideY - 1; // neighbor index (get convention) - m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ - nn = ijk - strideY + 1; // neighbor index (get convention) - m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ - nn = ijk - strideZ - 1; // neighbor index (get convention) - m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ - nn = ijk + strideZ + 1; // neighbor index (get convention) - m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ - nn = ijk + strideZ - 1; // neighbor index (get convention) - m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ - nn = ijk - strideZ + 1; // neighbor index (get convention) - m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ - nn = ijk - strideZ - strideY; // neighbor index (get convention) - m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ - nn = ijk + strideZ + strideY; // neighbor index (get convention) - m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ - nn = ijk + strideZ - strideY; // neighbor index (get convention) - m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ - nn = ijk - strideZ + strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 - //............Compute the Color Gradient................................... - nx = -(m1 - m2 + 0.5 * (m7 - m8 + m9 - m10 + m11 - m12 + m13 - m14)); - ny = -(m3 - m4 + 0.5 * (m7 - m8 - m9 + m10 + m15 - m16 + m17 - m18)); - nz = -(m5 - m6 + 0.5 * (m11 - m12 - m13 + m14 + m15 - m16 - m17 + m18)); - - //...........Normalize the Color Gradient................................. - C = sqrt(nx * nx + ny * ny + nz * nz); - double ColorMag = C; - if (C == 0.0) - ColorMag = 1.0; - nx = nx / ColorMag; - ny = ny / ColorMag; - nz = nz / ColorMag; + ijk = Map[n]; + + nn = ijk - 1; // neighbor index (get convention) + m1 = Phi[nn]; // get neighbor for phi - 1 + + nn = ijk + 1; // neighbor index (get convention) + m2 = Phi[nn]; // get neighbor for phi - 2 + + nn = ijk - strideY; // neighbor index (get convention) + m3 = Phi[nn]; // get neighbor for phi - 3 + + nn = ijk + strideY; // neighbor index (get convention) + m4 = Phi[nn]; // get neighbor for phi - 4 + + nn = ijk - strideZ; // neighbor index (get convention) + m5 = Phi[nn]; // get neighbor for phi - 5 + + nn = ijk + strideZ; // neighbor index (get convention) + m6 = Phi[nn]; // get neighbor for phi - 6 + + nn = ijk - strideY - 1; // neighbor index (get convention) + m7 = Phi[nn]; // get neighbor for phi - 7 + + nn = ijk + strideY + 1; // neighbor index (get convention) + m8 = Phi[nn]; // get neighbor for phi - 8 + + nn = ijk + strideY - 1; // neighbor index (get convention) + m9 = Phi[nn]; // get neighbor for phi - 9 + + nn = ijk - strideY + 1; // neighbor index (get convention) + m10 = Phi[nn]; // get neighbor for phi - 10 + + nn = ijk - strideZ - 1; // neighbor index (get convention) + m11 = Phi[nn]; // get neighbor for phi - 11 + + nn = ijk + strideZ + 1; // neighbor index (get convention) + m12 = Phi[nn]; // get neighbor for phi - 12 + + nn = ijk + strideZ - 1; // neighbor index (get convention) + m13 = Phi[nn]; // get neighbor for phi - 13 + + nn = ijk - strideZ + 1; // neighbor index (get convention) + m14 = Phi[nn]; // get neighbor for phi - 14 + + nn = ijk - strideZ - strideY; // neighbor index (get convention) + m15 = Phi[nn]; // get neighbor for phi - 15 + + nn = ijk + strideZ + strideY; // neighbor index (get convention) + m16 = Phi[nn]; // get neighbor for phi - 16 + + nn = ijk + strideZ - strideY; // neighbor index (get convention) + m17 = Phi[nn]; // get neighbor for phi - 17 + + nn = ijk - strideZ + strideY; // neighbor index (get convention) + m18 = Phi[nn]; // get neighbor for phi - 18 + + //............Compute the Color Gradient................................... + nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); + nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); + + //...........Normalize the Color Gradient................................. + C = sqrt(nx*nx+ny*ny+nz*nz); + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } + // q=0 fq = dist[n]; rho = fq; @@ -1989,7 +2151,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nx; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspx; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta; @@ -2005,7 +2167,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 2 // Cq = {0,1,0} - delta = beta * nA * nB * nAB * 0.1111111111111111 * ny; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspy; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta; @@ -2020,7 +2182,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 4 // Cq = {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nz; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspz; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta; @@ -2041,7 +2203,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // double Fx, double Fy, double Fz, int start, int finish, int Np){ extern "C" void ScaLBL_D3Q19_AAodd_Color( int *neighborList, int *Map, double *dist, double *Aq, double *Bq, - double *Den, double *Phi, double *Vel, double rhoA, double rhoB, + double *Den, double *Phi, unsigned int* NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np) { @@ -2061,6 +2223,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( double C, nx, ny, nz; //color gradient magnitude and direction double ux, uy, uz; double phi, tau, rho0, rlx_setA, rlx_setB; + double nspx, nspy, nspz; const double mrt_V1 = 0.05263157894736842; const double mrt_V2 = 0.012531328320802; @@ -2087,82 +2250,242 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = tauA + 0.5 * (1.0 - phi) * (tauB - tauA); + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); - // Get the 1D index based on regular data layout + // Get the 1D index based on regular data layout ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ + nn = ijk - 1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ + nn = ijk + 1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ + nn = ijk - strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ + nn = ijk + strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ + nn = ijk - strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ + nn = ijk + strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ + nn = ijk - strideY - 1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ + nn = ijk + strideY + 1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ + nn = ijk + strideY - 1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ + nn = ijk - strideY + 1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ + nn = ijk - strideZ - 1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ + nn = ijk + strideZ + 1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ + nn = ijk + strideZ - 1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ + nn = ijk - strideZ + 1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ + nn = ijk - strideZ - strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ + nn = ijk + strideZ + strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ + nn = ijk + strideZ - strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ + nn = ijk - strideZ + strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 + m18 = Phi[nn]; // get neighbor for phi - 18 + //............Compute the Color Gradient................................... - nx = -(m1 - m2 + 0.5 * (m7 - m8 + m9 - m10 + m11 - m12 + m13 - m14)); - ny = -(m3 - m4 + 0.5 * (m7 - m8 - m9 + m10 + m15 - m16 + m17 - m18)); - nz = -(m5 - m6 + 0.5 * (m11 - m12 - m13 + m14 + m15 - m16 - m17 + m18)); + nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); + nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); //...........Normalize the Color Gradient................................. - C = sqrt(nx * nx + ny * ny + nz * nz); - double ColorMag = C; - if (C == 0.0) - ColorMag = 1.0; - nx = nx / ColorMag; - ny = ny / ColorMag; - nz = nz / ColorMag; + C = sqrt(nx*nx+ny*ny+nz*nz); + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } // q=0 fq = dist[n]; @@ -2666,7 +2989,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nx; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspx; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta; @@ -2685,7 +3008,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // Cq = {0,1,0} - delta = beta * nA * nB * nAB * 0.1111111111111111 * ny; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspy; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta; @@ -2705,7 +3028,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // q = 4 // Cq = {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nz; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspz; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta; diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index 2a85dd40..dd37bec6 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -918,6 +918,51 @@ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np) { } } +extern "C" void ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int Np, double Fx, double Fy, double Fz) { + int n; + int N = Np; + // distributions + double f0, f1, f2, f3, f4, f5, f6, f7, f8, f9; + double f10, f11, f12, f13, f14, f15, f16, f17, f18; + double rho,vx, vy, vz; + + for (n = 0; n < N; n++) { + //........................................................................ + // Registers to store the distributions + //........................................................................ + f2 = dist[2 * N + n]; + f4 = dist[4 * N + n]; + f6 = dist[6 * N + n]; + f8 = dist[8 * N + n]; + f10 = dist[10 * N + n]; + f12 = dist[12 * N + n]; + f14 = dist[14 * N + n]; + f16 = dist[16 * N + n]; + f18 = dist[18 * N + n]; + //........................................................................ + f0 = dist[n]; + f1 = dist[N + n]; + f3 = dist[3 * N + n]; + f5 = dist[5 * N + n]; + f7 = dist[7 * N + n]; + f9 = dist[9 * N + n]; + f11 = dist[11 * N + n]; + f13 = dist[13 * N + n]; + f15 = dist[15 * N + n]; + f17 = dist[17 * N + n]; + //.................Compute the velocity................................... + rho = f0 + f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10 + f11 + f12 + f13 + f14 + f15 + f16 + f17 + f18; + vx = f1 - f2 + f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14; + vy = f3 - f4 + f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18; + vz = f5 - f6 + f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18; + //..................Write the velocity..................................... + vel[n] = vx/rho - Fx*0.5; + vel[N + n] = vy/rho - Fy*0.5; + vel[2 * N + n] = vz/rho - Fz*0.5; + //........................................................................ + } +} + extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *Pressure, int N) { for (int n = 0; n < N; n++) { //........................................................................ diff --git a/cuda/Color.cu b/cuda/Color.cu index c25e0f08..d4e8b90f 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -20,6 +20,7 @@ #define NBLOCKS 1024 #define NTHREADS 256 +#define f64_eps 1E-12 __global__ void dvc_ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz) { @@ -1266,7 +1267,7 @@ __global__ void dvc_ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int } -__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, +__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int *NeighborSolid, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; @@ -1281,6 +1282,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; + double nspx, nspy, nspz; // const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1311,69 +1313,67 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau=tauA + 0.5*(1.0-phi)*(tauB-tauA); + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); // Get the 1D index based on regular data layout ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ - nn = ijk-1; // neighbor index (get convention) - m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ - nn = ijk+1; // neighbor index (get convention) - m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ - nn = ijk-strideY; // neighbor index (get convention) - m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ - nn = ijk+strideY; // neighbor index (get convention) - m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ - nn = ijk-strideZ; // neighbor index (get convention) - m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ - nn = ijk+strideZ; // neighbor index (get convention) - m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ - nn = ijk-strideY-1; // neighbor index (get convention) - m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ - nn = ijk+strideY+1; // neighbor index (get convention) - m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ - nn = ijk+strideY-1; // neighbor index (get convention) - m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ - nn = ijk-strideY+1; // neighbor index (get convention) - m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ - nn = ijk-strideZ-1; // neighbor index (get convention) - m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ - nn = ijk+strideZ+1; // neighbor index (get convention) - m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ - nn = ijk+strideZ-1; // neighbor index (get convention) - m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ - nn = ijk-strideZ+1; // neighbor index (get convention) - m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ - nn = ijk-strideZ-strideY; // neighbor index (get convention) - m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ - nn = ijk+strideZ+strideY; // neighbor index (get convention) - m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ - nn = ijk+strideZ-strideY; // neighbor index (get convention) - m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ - nn = ijk-strideZ+strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 + + nn = ijk - 1; // neighbor index (get convention) + m1 = Phi[nn]; // get neighbor for phi - 1 + + nn = ijk + 1; // neighbor index (get convention) + m2 = Phi[nn]; // get neighbor for phi - 2 + + nn = ijk - strideY; // neighbor index (get convention) + m3 = Phi[nn]; // get neighbor for phi - 3 + + nn = ijk + strideY; // neighbor index (get convention) + m4 = Phi[nn]; // get neighbor for phi - 4 + + nn = ijk - strideZ; // neighbor index (get convention) + m5 = Phi[nn]; // get neighbor for phi - 5 + + nn = ijk + strideZ; // neighbor index (get convention) + m6 = Phi[nn]; // get neighbor for phi - 6 + + nn = ijk - strideY - 1; // neighbor index (get convention) + m7 = Phi[nn]; // get neighbor for phi - 7 + + nn = ijk + strideY + 1; // neighbor index (get convention) + m8 = Phi[nn]; // get neighbor for phi - 8 + + nn = ijk + strideY - 1; // neighbor index (get convention) + m9 = Phi[nn]; // get neighbor for phi - 9 + + nn = ijk - strideY + 1; // neighbor index (get convention) + m10 = Phi[nn]; // get neighbor for phi - 10 + + nn = ijk - strideZ - 1; // neighbor index (get convention) + m11 = Phi[nn]; // get neighbor for phi - 11 + + nn = ijk + strideZ + 1; // neighbor index (get convention) + m12 = Phi[nn]; // get neighbor for phi - 12 + + nn = ijk + strideZ - 1; // neighbor index (get convention) + m13 = Phi[nn]; // get neighbor for phi - 13 + + nn = ijk - strideZ + 1; // neighbor index (get convention) + m14 = Phi[nn]; // get neighbor for phi - 14 + + nn = ijk - strideZ - strideY; // neighbor index (get convention) + m15 = Phi[nn]; // get neighbor for phi - 15 + + nn = ijk + strideZ + strideY; // neighbor index (get convention) + m16 = Phi[nn]; // get neighbor for phi - 16 + + nn = ijk + strideZ - strideY; // neighbor index (get convention) + m17 = Phi[nn]; // get neighbor for phi - 17 + + nn = ijk - strideZ + strideY; // neighbor index (get convention) + m18 = Phi[nn]; // get neighbor for phi - 18 + //............Compute the Color Gradient................................... nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); @@ -1381,11 +1381,174 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //...........Normalize the Color Gradient................................. C = sqrt(nx*nx+ny*ny+nz*nz); - double ColorMag = C; - if (C==0.0) ColorMag=1.0; - nx = nx/ColorMag; - ny = ny/ColorMag; - nz = nz/ColorMag; + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } // q=0 fq = dist[n]; @@ -1660,7 +1823,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ - if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -19*alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); @@ -1802,7 +1964,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nx; + delta = beta*nA*nB*nAB*0.1111111111111111*nspx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -1817,7 +1979,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 2 // Cq = {0,1,0} - delta = beta*nA*nB*nAB*0.1111111111111111*ny; + delta = beta*nA*nB*nAB*0.1111111111111111*nspy; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -1831,7 +1993,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 4 // Cq = {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nz; + delta = beta*nA*nB*nAB*0.1111111111111111*nspz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -1850,7 +2012,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double *Phi, unsigned int *NeighborSolid, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ int n,nn,ijk,nread; @@ -1869,6 +2031,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; + double nspx, nspy, nspz; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1898,69 +2061,67 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau=tauA + 0.5*(1.0-phi)*(tauB-tauA); + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); // Get the 1D index based on regular data layout ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ - nn = ijk-1; // neighbor index (get convention) - m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ - nn = ijk+1; // neighbor index (get convention) - m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ - nn = ijk-strideY; // neighbor index (get convention) - m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ - nn = ijk+strideY; // neighbor index (get convention) - m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ - nn = ijk-strideZ; // neighbor index (get convention) - m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ - nn = ijk+strideZ; // neighbor index (get convention) - m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ - nn = ijk-strideY-1; // neighbor index (get convention) - m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ - nn = ijk+strideY+1; // neighbor index (get convention) - m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ - nn = ijk+strideY-1; // neighbor index (get convention) - m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ - nn = ijk-strideY+1; // neighbor index (get convention) - m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ - nn = ijk-strideZ-1; // neighbor index (get convention) - m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ - nn = ijk+strideZ+1; // neighbor index (get convention) - m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ - nn = ijk+strideZ-1; // neighbor index (get convention) - m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ - nn = ijk-strideZ+1; // neighbor index (get convention) - m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ - nn = ijk-strideZ-strideY; // neighbor index (get convention) - m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ - nn = ijk+strideZ+strideY; // neighbor index (get convention) - m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ - nn = ijk+strideZ-strideY; // neighbor index (get convention) - m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ - nn = ijk-strideZ+strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 + + nn = ijk - 1; // neighbor index (get convention) + m1 = Phi[nn]; // get neighbor for phi - 1 + + nn = ijk + 1; // neighbor index (get convention) + m2 = Phi[nn]; // get neighbor for phi - 2 + + nn = ijk - strideY; // neighbor index (get convention) + m3 = Phi[nn]; // get neighbor for phi - 3 + + nn = ijk + strideY; // neighbor index (get convention) + m4 = Phi[nn]; // get neighbor for phi - 4 + + nn = ijk - strideZ; // neighbor index (get convention) + m5 = Phi[nn]; // get neighbor for phi - 5 + + nn = ijk + strideZ; // neighbor index (get convention) + m6 = Phi[nn]; // get neighbor for phi - 6 + + nn = ijk - strideY - 1; // neighbor index (get convention) + m7 = Phi[nn]; // get neighbor for phi - 7 + + nn = ijk + strideY + 1; // neighbor index (get convention) + m8 = Phi[nn]; // get neighbor for phi - 8 + + nn = ijk + strideY - 1; // neighbor index (get convention) + m9 = Phi[nn]; // get neighbor for phi - 9 + + nn = ijk - strideY + 1; // neighbor index (get convention) + m10 = Phi[nn]; // get neighbor for phi - 10 + + nn = ijk - strideZ - 1; // neighbor index (get convention) + m11 = Phi[nn]; // get neighbor for phi - 11 + + nn = ijk + strideZ + 1; // neighbor index (get convention) + m12 = Phi[nn]; // get neighbor for phi - 12 + + nn = ijk + strideZ - 1; // neighbor index (get convention) + m13 = Phi[nn]; // get neighbor for phi - 13 + + nn = ijk - strideZ + 1; // neighbor index (get convention) + m14 = Phi[nn]; // get neighbor for phi - 14 + + nn = ijk - strideZ - strideY; // neighbor index (get convention) + m15 = Phi[nn]; // get neighbor for phi - 15 + + nn = ijk + strideZ + strideY; // neighbor index (get convention) + m16 = Phi[nn]; // get neighbor for phi - 16 + + nn = ijk + strideZ - strideY; // neighbor index (get convention) + m17 = Phi[nn]; // get neighbor for phi - 17 + + nn = ijk - strideZ + strideY; // neighbor index (get convention) + m18 = Phi[nn]; // get neighbor for phi - 18 + //............Compute the Color Gradient................................... nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); @@ -1968,11 +2129,174 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //...........Normalize the Color Gradient................................. C = sqrt(nx*nx+ny*ny+nz*nz); - double ColorMag = C; - if (C==0.0) ColorMag=1.0; - nx = nx/ColorMag; - ny = ny/ColorMag; - nz = nz/ColorMag; + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz ); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } // q=0 fq = dist[n]; @@ -2298,7 +2622,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ - if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -19*alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); @@ -2451,7 +2774,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nx; + delta = beta*nA*nB*nAB*0.1111111111111111*nspx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2469,7 +2792,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // Cq = {0,1,0} - delta = beta*nA*nB*nAB*0.1111111111111111*ny; + delta = beta*nA*nB*nAB*0.1111111111111111*nspy; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2488,7 +2811,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // q = 4 // Cq = {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nz; + delta = beta*nA*nB*nAB*0.1111111111111111*nspz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2877,7 +3200,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_ColorMomentum(int *neighborList, double * //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ - if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); @@ -4007,14 +4329,14 @@ extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A } // Pressure Boundary Conditions Functions -extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, +extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ cudaProfilerStart(); cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_Color, cudaFuncCachePreferL1); - dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB, + dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, NeighborSolid, Vel, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -4025,13 +4347,13 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do } extern "C" void ScaLBL_D3Q19_AAodd_Color(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ cudaProfilerStart(); cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_Color, cudaFuncCachePreferL1); - dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Vel, + dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, NeighborSolid, Vel, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); @@ -4164,4 +4486,3 @@ extern "C" void ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Sour dvc_ScaLBL_CopySlice_z<<>>(Phi,Nx,Ny,Nz,Source,Dest); } - diff --git a/cuda/D3Q19.cu b/cuda/D3Q19.cu index f895aff5..c7d8effc 100644 --- a/cuda/D3Q19.cu +++ b/cuda/D3Q19.cu @@ -1593,6 +1593,55 @@ __global__ void dvc_ScaLBL_D3Q19_Momentum(double *dist, double *vel, int N) } } +__global__ void dvc_ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int N, double Fx, double Fy, double Fz) +{ + int n; + // distributions + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + double rho,vx,vy,vz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(dist, vel, Np, Fx, Fy, Fz); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Velocity_2nd_order: %s \n",cudaGetErrorString(err)); + } +} + + extern "C" void ScaLBL_D3Q19_Pressure(double *fq, double *Pressure, int Np){ dvc_ScaLBL_D3Q19_Pressure<<< NBLOCKS,NTHREADS >>>(fq, Pressure, Np); } diff --git a/docs/source/examples/morphology/micp.rst b/docs/source/examples/morphology/micp.rst new file mode 100644 index 00000000..70ef0f36 --- /dev/null +++ b/docs/source/examples/morphology/micp.rst @@ -0,0 +1,29 @@ +************************************ +MICP Morphology Simulator +************************************ + +To expand and improve the current morphological tools in LBPM, new +morphology pre-processors and simulators were added. These changes +encompass the implementation of state-of-the-art morphology algorithms, +as well as an additional MICP simulator in which the invading non-wetting +phase is injected through all faces of the sample. + +.. code:: c + + Domain { + Filename = "crop_bt_101x112x88_uint8.raw" + N = 101, 112, 88 // domain size + n = 101, 112, 88 + nproc = 1, 1, 1 + ReadValues = 0, 1, 2 + WriteValues = 0, 1, 2 + voxel_length = 1 + BC = 0 + } + + FM { + Diameters = 1, 30, 1 // start, end, step + SaveImage = true + direction = "surround" // axis of intrusion + protocol = "micp" + } diff --git a/docs/source/userGuide/IO/fileformat.rst b/docs/source/userGuide/IO/fileformat.rst index c1da21cf..012a8088 100644 --- a/docs/source/userGuide/IO/fileformat.rst +++ b/docs/source/userGuide/IO/fileformat.rst @@ -2,12 +2,14 @@ I/O conventions for LBPM ======================== -There are three main kinds of output file that are supported by LBPM. +There are four main kinds of output file that are supported by LBPM. * CSV files -- space-delimited CSV files are used by the internal analysis framework * formatted binary files -- SILO and HDF5 formats are supported for visualization data + +* formatted binary files -- XML VTK format is also supported for visualization data * unformatted binary files -- ``.raw`` extension diff --git a/docs/source/userGuide/models/color/protocols/sw_steady.rst b/docs/source/userGuide/models/color/protocols/sw_steady.rst new file mode 100644 index 00000000..d63d066a --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/sw_steady.rst @@ -0,0 +1,72 @@ +====================================== +Color model -- Sw Steady +====================================== + +The water saturation steady state protocol is identical to the centrifuge protocol, with +the exception that the simulation explicity converges on the saturation state of fluid +A. + +That is, the simulation exits when + +.. math:: + :nowrap: + + $$ + \frac{\left | S_{w, i+1} - S_{w, i} \right |}{S_{w, i}} \le \epsilon + $$ + +averaged over the ``analysis_interval`` where :math:`S_{w,i}` is the saturation of fluid +A at step :math:`i` and :math:`\epsilon` is an allowed convergence threshold, or when +the ``timestepMax`` is reached, whichever occurs first. + +By default, :math:`\epsilon` is set to zero (i.e., the simulation continues through the +``timestepMax``), but can be set via ``tolerance`` within the ``Analysis`` section of +the input file database as shown below. + + +.. code-block:: c + + Color { + protocol = "sw_steady" + timestepMax = 1000000 // maximum timtestep + alpha = 0.005 // controls interfacial tension + rhoA = 1.0 // controls the density of fluid A + rhoB = 1.0 // controls the density of fluid B + tauA = 0.7 // controls the viscosity of fluid A + tauB = 0.7 // controls the viscosity of fluid B + F = 0, 0, -1.0e-5 // body force + din = 1.0 // inlet density (controls pressure) + dout = 1.0 // outlet density (controls pressure) + WettingConvention = "SCAL" // convention for sign of wetting affinity + ComponentLabels = 0, -1, -2 // image labels for solid voxels + ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label + Restart = false + } + Domain { + Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw" + ReadType = "8bit" // data type + N = 900, 900, 1600 // size of original image + nproc = 2, 2, 2 // process grid + n = 200, 200, 200 // sub-domain size + offset = 300, 300, 300 // offset to read sub-domain + voxel_length = 1.66 // voxel length (in microns) + ReadValues = -2, -1, 0, 1, 2 // labels within the original image + WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM + BC = 3 // boundary condition type (0 for periodic) + } + Analysis { + analysis_interval = 1000 // logging interval for timelog.csv + subphase_analysis_interval = 5000 // loggging interval for subphase.csv + visualization_interval = 100000 // interval to write visualization files + N_threads = 4 // number of analysis threads (GPU version only) + restart_interval = 1000000 // interval to write restart file + restart_file = "Restart" // base name of restart file + tolerance = 1e-9 // Sw convergence tolerance + } + Visualization { + write_silo = true // write SILO databases with assigned variables + save_8bit_raw = true // write labeled 8-bit binary files with phase assignments + save_phase_field = true // save phase field within SILO database + save_pressure = false // save pressure field within SILO database + save_velocity = false // save velocity field within SILO database + } diff --git a/docs/source/userGuide/visualization/vtk.rst b/docs/source/userGuide/visualization/vtk.rst new file mode 100644 index 00000000..69a519f9 --- /dev/null +++ b/docs/source/userGuide/visualization/vtk.rst @@ -0,0 +1,19 @@ +====================================== +VTK XML output format +====================================== + +Due to several limitations in the support for VisIt files in ParaView, we added support for VTK XML (.vti) files, which are natively supported by ParaView. + +This output format can be enabled by setting format = "vtk" in the Visualization section of the database file. The fields to be written to output can be chosen using the same parameters as in the VisIt output. + +.. code:: c + + Visualization { + format = "vtk" + save_phase_field = true + save_pressure = false + } + +The interval for writing output files is defined by the ``visualization_interval`` parameter in the ``Analysis`` section. + +LBPM also generates and updates a ``LBM.pvd`` file, which allows ParaView to open all .vti files generated during the current simulation. diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 88caa2ea..aedc4795 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -32,7 +32,7 @@ ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, tauB(0), rhoA(0), rhoB(0), alpha(0), beta(0), Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), inletA(0), inletB(0), outletA(0), outletB(0), Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), - BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), + BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), NeighborSolid(nullptr), NeighborList(nullptr), dvcMap(nullptr), fq(nullptr), Aq(nullptr), Bq(nullptr), Den(nullptr), Phi(nullptr), ColorGrad(nullptr), Velocity(nullptr), Pressure(nullptr), comm(COMM) { @@ -111,6 +111,9 @@ void ScaLBL_ColorModel::ReadParams(string filename) { if (color_db->keyExists("flux")) { flux = color_db->getScalar("flux"); } + if (color_db->keyExists("timestep")) { + timestep = color_db->getScalar("timestep"); + } inletA = 1.f; inletB = 0.f; outletA = 0.f; @@ -180,7 +183,15 @@ void ScaLBL_ColorModel::ReadParams(string filename) { "pressure boundary condition \n"); } domain_db->putScalar("BC", BoundaryCondition); - } else if (protocol == "core flooding") { + } else if (protocol == "sw_steady") { + if (BoundaryCondition != 3) { + BoundaryCondition = 3; + if (rank == 0) + printf("WARNING: protocol (sw_steady) supports only constant " + "pressure boundary condition \n"); + } + domain_db->putScalar("BC", BoundaryCondition); + } else if (protocol == "core flooding") { if (rank == 0) printf("Using core flooding protocol \n"); if (BoundaryCondition != 4) { @@ -448,6 +459,8 @@ void ScaLBL_ColorModel::Create() { ScaLBL_AllocateDeviceMemory((void **)&Pressure, sizeof(double) * Np); ScaLBL_AllocateDeviceMemory((void **)&Velocity, 3 * sizeof(double) * Np); ScaLBL_AllocateDeviceMemory((void **)&ColorGrad, 3 * sizeof(double) * Np); + ScaLBL_AllocateDeviceMemory((void **)&NeighborSolid, sizeof(unsigned int) * Np); + //........................................................................... // Update GPU data structures if (rank == 0) @@ -503,6 +516,43 @@ void ScaLBL_ColorModel::Create() { if (rank == 0) printf("Model created \n"); delete[] PhaseLabel; + + unsigned int *TmpSolid = new unsigned int[Np]; + + for (int k = 1; k < Nz - 1; k++) { + for (int j = 1; j < Ny - 1; j++) { + for (int i = 1; i < Nx - 1; i++) { + int idx = Map(i, j, k); + unsigned int data = 0; + if (!(idx < 0)) { + if ((Map(i - 1 , j, k) == -1) ) data |= (1u << 1); + if ((Map(i + 1 , j, k) == -1) ) data |= (1u << 2); + if ((Map(i , j - 1, k) == -1) ) data |= (1u << 3); + if ((Map(i , j + 1, k) == -1) ) data |= (1u << 4); + if ((Map(i , j, k - 1) == -1) ) data |= (1u << 5); + if ((Map(i , j, k + 1) == -1) ) data |= (1u << 6); + if ((Map(i - 1 , j - 1, k) == -1) ) data |= (1u << 7); + if ((Map(i + 1 , j + 1, k) == -1) ) data |= (1u << 8); + if ((Map(i - 1 , j + 1, k) == -1) ) data |= (1u << 9); + if ((Map(i + 1, j - 1, k) == -1) ) data |= (1u << 10); + if ((Map(i - 1 , j, k - 1) == -1) ) data |= (1u << 11); + if ((Map(i + 1 , j, k + 1) == -1) ) data |= (1u << 12); + if ((Map(i - 1 , j, k + 1) == -1) ) data |= (1u << 13); + if ((Map(i + 1 , j, k - 1) == -1) ) data |= (1u << 14); + if ((Map(i , j - 1, k - 1) == -1) ) data |= (1u << 15); + if ((Map(i , j + 1, k + 1) == -1) ) data |= (1u << 16); + if ((Map(i , j - 1, k + 1) == -1) ) data |= (1u << 17); + if ((Map(i , j + 1, k - 1) == -1) ) data |= (1u << 18); + TmpSolid[idx] = data; + } + } + } + } + + ScaLBL_CopyToDevice(NeighborSolid, TmpSolid, sizeof(unsigned int) * Np); + ScaLBL_Comm->Barrier(); + delete[] TmpSolid; + } /******************************************************** @@ -630,7 +680,6 @@ double ScaLBL_ColorModel::Run(int returntime) { bool RESCALE_FORCE = false; bool SET_CAPILLARY_NUMBER = false; bool TRIGGER_FORCE_RESCALE = false; - double tolerance = 0.01; auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); auto current_db = db->cloneDatabase(); auto flow_db = db->getDatabase("FlowAdaptor"); @@ -656,9 +705,8 @@ double ScaLBL_ColorModel::Run(int returntime) { color_db->getScalar("rescale_force_after_timestep"); RESCALE_FORCE = true; } - if (analysis_db->keyExists("tolerance")) { - tolerance = analysis_db->getScalar("tolerance"); - } + double tolerance = analysis_db->getWithDefault("tolerance", 1e-5); + int analysis_interval = analysis_db->getWithDefault("analysis_interval", 1000); runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map); @@ -690,7 +738,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); @@ -709,7 +757,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); @@ -735,7 +783,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -754,7 +802,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); @@ -764,7 +812,7 @@ double ScaLBL_ColorModel::Run(int returntime) { Den); // allow initial ramp-up to get closer to steady state CURRENT_TIMESTEP += 2; - if (CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { + if (CURRENT_TIMESTEP % analysis_interval == 0 && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { analysis.finish(); double volB = Averages->gwb.V; @@ -801,7 +849,7 @@ double ScaLBL_ColorModel::Run(int returntime) { fabs(muA * flow_rate_A + muB * flow_rate_B) / (5.796 * alpha); bool isSteady = false; - if ((fabs((Ca - Ca_previous) / Ca) < tolerance && + if ((fabs((Ca - Ca_previous) / analysis_interval / Ca) < tolerance && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS)) isSteady = true; if (CURRENT_TIMESTEP >= MAX_STEADY_TIMESTEPS) @@ -1071,7 +1119,12 @@ double ScaLBL_ColorModel::Run(int returntime) { printf("Ca = %f, (previous = %f) \n", Ca, Ca_previous); } } + + break; // steady-state achieved, exit. } + + // save for convergence checks + Ca_previous = Ca; } } analysis.finish(); @@ -1100,11 +1153,9 @@ double ScaLBL_ColorModel::Run(int returntime) { void ScaLBL_ColorModel::Run() { int nprocs = nprocx * nprocy * nprocz; const RankInfoStruct rank_info(rank, nprocx, nprocy, nprocz); - int analysis_interval = - 1000; // number of timesteps in between in situ analysis - if (analysis_db->keyExists("analysis_interval")) { - analysis_interval = analysis_db->getScalar("analysis_interval"); - } + + int analysis_interval = analysis_db->getWithDefault("analysis_interval", 1000); + double tolerance = analysis_db->getWithDefault("tolerance", 0.0); //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); @@ -1116,7 +1167,11 @@ void ScaLBL_ColorModel::Run() { Map); //analysis.createThreads( analysis_method, 4 ); auto t1 = std::chrono::system_clock::now(); - while (timestep < timestepMax) { + + double delta_sw = 1.0; + double sw_prev = -1.0; + + while (timestep < timestepMax && delta_sw > tolerance) { PROFILE_START("Update"); // *************ODD TIMESTEP************* @@ -1142,7 +1197,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); @@ -1161,7 +1216,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); @@ -1187,7 +1242,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -1206,20 +1261,30 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); //************************************************************************ PROFILE_STOP("Update"); - if (rank == 0 && timestep % analysis_interval == 0 && - BoundaryCondition == 4) { - printf("%i %f \n", timestep, din); - } // Run the analysis analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den); + + if (timestep % analysis_interval == 0){ + analysis.finish(); + + double volA = Averages->gnb.V / Dm->Volume; + double volB = Averages->gwb.V / Dm->Volume; + double sw = volB / (volA + volB); + + delta_sw = fabs(sw - sw_prev) / analysis_interval / sw; + if (rank == 0) + printf("t: %d sw: %0.5e dSw/dt: %.5e\n", timestep, sw, delta_sw); + + sw_prev = sw; + } } analysis.finish(); PROFILE_STOP("Loop"); diff --git a/models/ColorModel.h b/models/ColorModel.h index 888f4fc2..b8a1dea5 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -138,6 +138,7 @@ class ScaLBL_ColorModel { IntArray Map; signed char *id; + unsigned int *NeighborSolid; int *NeighborList; int *dvcMap; double *fq, *Aq, *Bq; diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index ae3f530e..674132f1 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -263,7 +263,7 @@ void ScaLBL_MRTModel::Run() { if (WriteHeader) { log_file = fopen("Permeability.csv", "a+"); - fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz absperm\n"); + fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz absperm(mDa) absperm*(mDa)\n"); fclose(log_file); } } @@ -328,7 +328,7 @@ void ScaLBL_MRTModel::Run() { //************************************************************************/ if (timestep % ANALYSIS_INTERVAL == 0) { - ScaLBL_D3Q19_Momentum(fq, Velocity, Np); + ScaLBL_D3Q19_Momentum_2nd_order(fq, Velocity, Np, Fx, Fy, Fz); ScaLBL_DeviceBarrier(); comm.barrier(); ScaLBL_Comm->RegularLayout(Map, &Velocity[0], Velocity_x); @@ -391,8 +391,7 @@ void ScaLBL_MRTModel::Run() { Xs = Dm->Comm.sumReduce(Xs); double h = Dm->voxel_length; - double absperm = - h * h * mu * Mask->Porosity() * Mask->Porosity() * flow_rate / force_mag; + double absperm = h * h * mu * Mask->Porosity() * flow_rate / force_mag; absperm *= 1013.0; // Convert to mDarcy if (rank == 0) { @@ -400,9 +399,9 @@ void ScaLBL_MRTModel::Run() { FILE *log_file = fopen("Permeability.csv", "a"); fprintf(log_file, "%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g " - "%.8g %.8g\n", + "%.8g %.8g %.8g\n", timestep, Fx, Fy, Fz, mu, h * h * h * Vs, h * h * As, - h * Hs, Xs, vax, vay, vaz, absperm); + h * Hs, Xs, vax, vay, vaz, absperm, absperm * Mask->Porosity()); fclose(log_file); } } @@ -484,7 +483,7 @@ void ScaLBL_MRTModel::VelocityField() { auto VzVar = std::make_shared(); auto SignDistVar = std::make_shared(); - IO::initialize("", format, "false"); + IO::initialize("", format, false); // Create the MeshDataStruct visData.resize(1); visData[0].meshName = "domain"; @@ -527,7 +526,7 @@ void ScaLBL_MRTModel::VelocityField() { fillData.copy(Velocity_x, VelxData); fillData.copy(Velocity_y, VelyData); fillData.copy(Velocity_z, VelzData); - + IO::writeData(timestep, visData, Dm->Comm); } } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ceafff4..f8a31b60 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -48,7 +48,7 @@ ADD_LBPM_EXECUTABLE( TestIonModel ) ADD_LBPM_EXECUTABLE( TestNernstPlanck ) ADD_LBPM_EXECUTABLE( TestPNP_Stokes ) ADD_LBPM_EXECUTABLE( TestMixedGrad ) - +ADD_LBPM_EXECUTABLE( fmm_micp_drainage ) CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY ) @@ -84,7 +84,7 @@ ADD_LBPM_TEST_1_2_4( TestBlobIdentify ) #ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 ) #ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 ) ADD_LBPM_TEST_PARALLEL( TestSegDist 8 ) -ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 ) +ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 test.db) ADD_LBPM_TEST_1_2_4( testCommunication ) ADD_LBPM_TEST( TestWriter ) ADD_LBPM_TEST( TestDatabase ) diff --git a/tests/DataAggregator.cpp b/tests/DataAggregator.cpp index 137ebd61..466b25a9 100644 --- a/tests/DataAggregator.cpp +++ b/tests/DataAggregator.cpp @@ -5,6 +5,7 @@ #include #include #include +#include using namespace std; diff --git a/tests/fmm_micp_drainage.cpp b/tests/fmm_micp_drainage.cpp new file mode 100644 index 00000000..ff37b057 --- /dev/null +++ b/tests/fmm_micp_drainage.cpp @@ -0,0 +1,727 @@ +/* +<<<<<<< HEAD + Copyright 2026 + Diogo Nardelli Siebert, Universidade Federal de Santa Catarina + Bernardo Gehlen, Universidade Federal de Santa Catarina + Alexandre Miers Zabot, Universidade Federal de Santa Catarina + + 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 :contentReference[oaicite:0]{index=0}, 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 . +*/ + +/* + ---------------------------------------------------------------------- + Methodological Notes + + This file implements the method described in: + + A. M. Zabot et al., + "A Unified Algorithm for the Young–Laplace Method Applied to Porous Media," + Brazilian Journal of Physics, vol. 54, no. 3, 2024, p. 63. + + The distance transform is computed using the method from: + + P. F. Felzenszwalb and D. P. Huttenlocher, + "Distance transforms of sampled functions," + Theory of Computing, vol. 8, no. 1, pp. 415–428, 2012. + + Connected component labeling is based on: + + L. He, Y. Chao, and K. Suzuki, + "A run-based two-scan labeling algorithm," + IEEE Transactions on Image Processing, vol. 17, no. 5, pp. 749–756, 2008. + ---------------------------------------------------------------------- +*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include "../common/Array.h" +#include "../common/Domain.h" +#include "../common/UtilityMacros.h" + +#define SOLID ((unsigned char)0) +#define DISPLACED ((unsigned char)1) +#define INJECTED ((unsigned char)2) + +#define BACKGROUND 0 +#define FOREGROUND 1 + +using namespace std; + +void merge(const int &u, const int &v, vector &next, vector &tail, + vector &rtable) { + for (int i = v; i != -1;) { + rtable[i] = u; + i = next[i]; + } + next[tail[u]] = v; + tail[u] = tail[v]; +} + +void resolve(const int &x, const int &y, vector &next, vector &tail, + vector &rtable) { + const int u = rtable[x]; + const int v = rtable[y]; + if (u < v) + merge(u, v, next, tail, rtable); + else if (v < u) + merge(v, u, next, tail, rtable); +} + +void component_labeling(IntArray &IMG, const int &F, const int &B) { + + size_t maxNumberOfLabels = (IMG.length() + 1) / 2; + + std::vector next(maxNumberOfLabels); + std::vector tail(maxNumberOfLabels); + std::vector rtable(maxNumberOfLabels); + + const int nx = IMG.size(0); + const int ny = IMG.size(1); + const int nz = IMG.size(2); + + int lx = 0, nl = 1; + vector uniq_labels(3); + int nuniq; + + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + + if (IMG(x, y, z) == F) { + const int lq = (x > 0) ? IMG(x - 1, y, z) : B; + const int lp = (y > 0) ? IMG(x, y - 1, z) : B; + const int lz = (z > 0) ? IMG(x, y, z - 1) : B; + + nuniq = 0; + if (lp != B) { + uniq_labels[nuniq] = lp; + nuniq++; + } + if (lq != B && lq != lp) { + uniq_labels[nuniq] = lq; + nuniq++; + } + if (lz != B && lz != lp && lz != lq) { + uniq_labels[nuniq] = lz; + nuniq++; + } + + // Deals with non unique labels in neighboaring points by unifying them + switch (nuniq) { + + case 0: + nl++; + lx = nl; + + rtable[nl] = nl; + next[nl] = -1; + tail[nl] = nl; + break; + + case 1: + lx = uniq_labels[0]; + break; + + case 2: + + resolve(uniq_labels[0], uniq_labels[1], next, tail, + rtable); + + lx = min(uniq_labels[0], uniq_labels[1]); + break; + + case 3: + + resolve(uniq_labels[0], uniq_labels[1], next, tail, + rtable); + resolve(uniq_labels[0], uniq_labels[2], next, tail, + rtable); + resolve(uniq_labels[1], uniq_labels[2], next, tail, + rtable); + + lx = min(uniq_labels[0], + min(uniq_labels[1], uniq_labels[2])); + break; + } + + IMG(x, y, z) = lx; + } + } + } + } + + int *img = IMG.data(); + + for (size_t n = 0; n < IMG.length(); n++) { + if (*img != B) + *img = rtable[*img]; + img++; + } +} + +// Compute the intersection between two lower-envelope parabolas. +static inline float intersection(int q, int vk, float fq, float fvk) { + float qq = (float)q * (float)q; + float vv = (float)vk * (float)vk; + return ((fq + qq) - (fvk + vv)) / (2.0 * ((float)q - (float)vk)); +} + +// Compute the exact squared Euclidean distance transform for a 1D line. +void edt_1d(const int *f, int *g, int n) { + int *v = (int *)malloc((size_t)n * sizeof(int)); + float *z = (float *)malloc((size_t)(n + 1) * sizeof(float)); + + int k = 0; + v[0] = 0; + z[0] = -INFINITY; + z[1] = INFINITY; + + for (int q = 1; q < n; q++) { + float s; + + while (1) { + int vk = v[k]; + s = intersection(q, vk, f[q], f[vk]); + + if (s <= z[k]) { + k--; + if (k < 0) { + k = 0; + break; + } + } else { + break; + } + } + + if (k == 0) { + int vk = v[k]; + s = intersection(q, vk, f[q], f[vk]); + + if (s <= z[k]) { + v[0] = q; + z[0] = -INFINITY; + z[1] = INFINITY; + continue; + } + } + + k++; + v[k] = q; + z[k] = s; + z[k + 1] = INFINITY; + } + + k = 0; + for (int x = 0; x < n; x++) { + while (z[k + 1] < (float)x) { + k++; + } + + int vk = v[k]; + int dx = x - vk; + g[x] = dx * dx + f[vk]; + } + + free(v); + free(z); +} + +// Copy one strided line, transform it in 1D, and write it back. +void process_line(int *edt2, const size_t first, const size_t stride, int n) { + + IntArray f(n); + IntArray g(n); + + size_t pos = first; + for (int i = 0; i < n; i++) { + f(i) = edt2[pos]; + pos += stride; + } + + edt_1d(f.data(), g.data(), n); + + pos = first; + for (int i = 0; i < n; i++) { + edt2[pos] = g(i); + pos += stride; + } +} + +template +void edt_3d(unsigned char target, Array &image, IntArray &distance2) { + const int nx = image.size(0); + const int ny = image.size(1); + const int nz = image.size(2); + + size_t nvox = image.length(); + int BIG = static_cast(nx * nx + ny * ny + nz * nz) + 1; + + TYPE *img = image.data(); + int *edt2 = distance2.data(); + + // Initialize target voxels with zero distance and all others with a large value. + for (int i = 0; i < (int)nvox; i++) { + edt2[i] = (img[i] == target) ? 0 : BIG; + } + + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + process_line(edt2, nx * (y + z * ny), 1, nx); + } + } + + for (int z = 0; z < nz; z++) { + for (int x = 0; x < nx; x++) { + process_line(edt2, z * nx * ny + x, nx, ny); + } + } + + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + process_line(edt2, y * nx + x, nx * ny, nz); + } + } +} + +// Validate database options and report the accepted values on failure. +void checkOption(std::string a, std::vector s, std::string keyName) { + std::string message = "Error: Invalid option '" + a + "' for " + keyName + + ". Valid options are: "; + for (string b : s) { + if (a == b) + return; + message += "'" + b + "', "; + } + message.pop_back(); + ERROR(message); +} + +template +void setRegion(Array &A, TYPE value, int x0, int x1, int y0, int y1, + int z0, int z1) { + for (int z = z0; z < z1; z++) + for (int y = y0; y < y1; y++) + for (int x = x0; x < x1; x++) { + A(x, y, z) = value; + } +} + +class Full_Morphology { +public: + Full_Morphology(int, char *[]); + int calc(const int &); + +public: + std::vector r_ini = {0, 0, 0}; + std::vector r_end = {0, 0, 0}; + + std::vector flowAxis = { + false, false, false}; // Injection occurs along the selected axis + bool flowPos = false; // Sense of invasion + + int ny, nx, nz; // Work dimensions (with reservoirs, if any) + int dimy, dimx, dimz; // Original dimensions + int NP; // Porous pixels + int inletPos, outletPos; // Reservoir regions + + int rChamberI[3]; + int rChamberO[3]; + + double resolution; + + vector diameter; // Invasion diameters + + bool compressible = false; // Compressibility, enabled for MICP + + bool allFaces; // Add reservoirs around all active faces for MICP + bool saveImg; + + Array originalState; // Original image, used for reference + Array currentState; // Work image, modified during the simulation + Array finalMap; + IntArray originalEDT; + Array trapped; +}; + +Full_Morphology::Full_Morphology(int argc, char *argv[]) { + + string filename; + + filename = argv[1]; + + auto db = std::make_shared(filename); + + auto domain_db = db->getDatabase("Domain"); + auto fm_db = db->getDatabase("FM"); + + auto size = domain_db->getVector("N"); + nx = size[0]; + ny = size[1]; + nz = size[2]; + + finalMap.resize(size[0], size[1], size[2]); + finalMap.fill(-1); + + auto ReadValues = domain_db->getVector("ReadValues"); + auto WriteValues = domain_db->getVector("WriteValues"); + + resolution = domain_db->getScalar("voxel_length"); + + auto READFILE = domain_db->getScalar("Filename"); + const string mmfile(READFILE); + + saveImg = fm_db->getScalar("SaveImage"); + + auto protocol = fm_db->getScalar("protocol"); + checkOption(protocol, {"micp", "drainage"}, "protocol"); + + if (protocol == "micp") { + compressible = true; + allFaces = true; + flowPos = true; + if (size[2] > 1) + flowAxis[2] = true; + else if (size[1] > 1) + flowAxis[1] = true; + else if (size[0] > 1) + flowAxis[0] = true; + } else if (protocol == "drainage") { + allFaces = false; + auto direction = fm_db->getScalar("direction"); + checkOption(direction, {"+x", "-x", "+y", "-y", "+z", "-z"}, + "direction"); + flowAxis[0] = (direction[1] == 'x'); + flowAxis[1] = (direction[1] == 'y'); + flowAxis[2] = (direction[1] == 'z'); + flowPos = (direction[0] == '+'); + compressible = fm_db->getWithDefault("compressible", false); + } + + auto diameterRange = fm_db->getVector("Diameters"); + int numberOfDiameters = + (diameterRange[1] - diameterRange[0]) / diameterRange[2] + 1; + + if (numberOfDiameters <= 0) { + ERROR("Error: It was impossible to create diameters array. "); + } + + diameter.resize(numberOfDiameters); + + for (int i = 0; i < numberOfDiameters; i++) { + diameter[i] = diameterRange[1] - i * diameterRange[2]; + } + + r_end = size; + dimx = nx; + dimy = ny; + dimz = nz; + + // Add additional layers for input/output reservoirs. + for (int i = 0; i < 3; i++) { + if (((flowAxis[i]) && (!allFaces)) || ((size[i] > 1) && (allFaces))) { + size[i] += 2; + r_ini[i] = 1; + r_end[i] = size[i] - 1; + } + } + + nx = size[0]; + ny = size[1]; + nz = size[2]; + + originalState.resize(nx, ny, nz); + currentState.resize(originalState.size()); + originalEDT.resize(originalState.size()); + trapped.resize(originalState.size()); + + currentState.fill(INJECTED); + trapped.fill(false); + + // For directional injection, set the outlet face to DISPLACED fluid. + for (int i = 0; i < 3; i++) { + int rMin[3] = {0, 0, 0}; + int rMax[3] = {nx, ny, nz}; + + if (!allFaces) { + if (flowAxis[i]) { + inletPos = flowPos ? 0 : size[i] - 1; + outletPos = flowPos ? size[i] - 1 : 0; + rMin[i] = outletPos; + rMax[i] = rMin[i] + 1; + setRegion(currentState, DISPLACED, rMin[0], rMax[0], rMin[1], + rMax[1], rMin[2], rMax[2]); + } + rChamberI[i] = flowAxis[i] ? inletPos : rMax[i] / 2; + rChamberO[i] = flowAxis[i] ? outletPos : rMax[i] / 2; + } + + else + rChamberI[i] = 0; + } + + int mapValue[255] = {-1}; + for (size_t idx = 0; idx < ReadValues.size(); idx++) { + + if ((ReadValues[idx] < 0) || (ReadValues[idx] > 255)) { + ERROR("Only values between 0 - 255 can be used as labels in " + "ReadValues"); + cout << ReadValues[idx] << endl; + } + if ((WriteValues[idx] < 0) || (WriteValues[idx] > 2)) { + ERROR("Only values between 0 (SOLID), 1 and 2 (INJECT/DISPLACED " + "FLUIDS) can be used as labels in WriteValues"); + } + mapValue[ReadValues[idx]] = (int)WriteValues[idx]; + } + + FILE *rawFile = fopen(mmfile.c_str(), "r"); + if (rawFile == NULL) { + ERROR("Error openning the file " + mmfile); + } + + long SEEK_BEGIN = ftell(rawFile); + long expectedSize = (long)(r_end[2] - r_ini[2]) * + (long)(r_end[1] - r_ini[1]) * + (long)(r_end[0] - r_ini[0]); + + fseek(rawFile, 0, SEEK_END); + + if (ftell(rawFile) != expectedSize) { + ERROR("File '" + mmfile + "' size is different from the expected (" + + to_string(expectedSize) + " bytes)."); + } + + fseek(rawFile, 0, + SEEK_BEGIN); // Move to the beginning of the file before reading. + + unsigned char readValue; + + NP = 0; + for (int z = r_ini[2]; z < r_end[2]; z++) { + for (int y = r_ini[1]; y < r_end[1]; y++) { + for (int x = r_ini[0]; x < r_end[0]; x++) { + + fread(&readValue, sizeof(unsigned char), 1, rawFile); + if (mapValue[readValue] == -1) { + ERROR(std::string("Not specified value in '" + filename + + "' at (" + to_string(x) + ", " + + to_string(y) + ", " + to_string(z) + + ").")); + } + + currentState(x, y, z) = (unsigned char)mapValue[readValue]; + if (currentState(x, y, z) != SOLID) + NP++; + } + } + } + + fclose(rawFile); + originalState = currentState; + + // Calculate the distance transform of the original solid phase. + edt_3d(SOLID, originalState, originalEDT); + + // Create the output CSV file and header when needed. + bool WriteHeader = false; + FILE *log_file = fopen("injection_output.csv", "r"); + if (log_file != NULL) + fclose(log_file); + else + WriteHeader = true; + + if (WriteHeader) { + log_file = fopen("injection_output.csv", "a+"); + fprintf(log_file, "step diameter_px diameter_um num_px_in frac_in " + "num_px_out frac_out\n"); + fclose(log_file); + } +} + +int Full_Morphology::calc(const int &step) { + + const int D = diameter[step]; + const double D24 = D * D / 4.0; + + IntArray auxMatrix(nx, ny, nz); + + // Mark the candidate invaded region as the region where the center of a + // D diameter sphere can be placed (erosion of the solid region) united with the reservoir. + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + auxMatrix(x, y, z) = (originalEDT(x, y, z) >= D24) || (originalState(x, y, z) == INJECTED) ? FOREGROUND : BACKGROUND; + } + } + } + + component_labeling(auxMatrix, FOREGROUND, BACKGROUND); + + // Extract the label associated with the injection layer. + int chamber_label = auxMatrix(rChamberI[0], rChamberI[1], rChamberI[2]); + + // Filter only the voxels connected to the input where the center of a a sphere with the specified diameter + // can be placed. + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + auxMatrix(x, y, z) = (originalEDT(x, y, z) >= D24 && auxMatrix(x,y,z) == chamber_label) ? FOREGROUND : BACKGROUND; + } + } + } + + // Compute the distance from the filtered region + edt_3d(FOREGROUND, auxMatrix, auxMatrix); + + // Perform the dilation of the filtered region and performs the union of this result + // with the previous step stored at the currentState + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + if (auxMatrix(x, y, z) < D24) { + currentState(x, y, z) = INJECTED; + } + auxMatrix(x, y, z) = (currentState(x, y, z) == INJECTED) + ? FOREGROUND + : BACKGROUND; + } + } + } + + + component_labeling(auxMatrix, FOREGROUND, BACKGROUND); + chamber_label = auxMatrix(rChamberI[0], rChamberI[1], rChamberI[2]); + + // Readds the DISPLACED fluid output layer + if (!allFaces) { + + if (flowAxis[0]) { + setRegion(currentState, DISPLACED, outletPos, outletPos + 1, 0, ny, + 0, nz); + } else if (flowAxis[1]) { + setRegion(currentState, DISPLACED, 0, nx, outletPos, outletPos + 1, + 0, nz); + } else if (flowAxis[2]) { + setRegion(currentState, DISPLACED, 0, nx, 0, ny, outletPos, + outletPos + 1); + } + } + + // For incompressible flow, disconnected displaced regions remain trapped + // in the final state. + if (!compressible) { + + for (int x = 0; x < nx; x++) { + for (int y = 0; y < ny; y++) { + for (int z = 0; z < nz; z++) { + + if (trapped(x, y, z)) + currentState(x, y, z) = DISPLACED; + + auxMatrix(x, y, z) = (currentState(x, y, z) == DISPLACED) + ? FOREGROUND + : BACKGROUND; + } + } + } + + component_labeling(auxMatrix, FOREGROUND, BACKGROUND); + chamber_label = auxMatrix(rChamberO[0], rChamberO[1], rChamberO[2]); + + for (int x = 0; x < nx; x++) { + for (int y = 0; y < ny; y++) { + for (int z = 0; z < nz; z++) { + if (currentState(x, y, z) == DISPLACED && + auxMatrix(x, y, z) != chamber_label) + trapped(x, y, z) = true; + } + } + } + } + + int injectedVolume = 0, displacedVolume = 0; + for (int z = r_ini[2]; z < r_end[2]; z++) { + for (int y = r_ini[1]; y < r_end[1]; y++) { + for (int x = r_ini[0]; x < r_end[0]; x++) { + + if (currentState(x, y, z) == INJECTED) + injectedVolume++; + else if (currentState(x, y, z) == DISPLACED) + displacedVolume++; + + int16_t *value = + &finalMap(x - r_ini[0], y - r_ini[1], z - r_ini[2]); + if (*value == -1 && currentState(x, y, z) == INJECTED) + *value = (int16_t)D; + if (*value == -1 && originalState(x, y, z) == SOLID) + *value = 0; + } + } + } + + if (saveImg && (D == diameter.back())) { + FILE *FRAW; + + FRAW = fopen("invasion_diameters.raw", "wb"); + fwrite(finalMap.data(), sizeof(int16_t), finalMap.length(), FRAW); + fclose(FRAW); + + FILE *FMHD = fopen("invasion_diameters.mhd", "w"); + + fprintf(FMHD, "ObjectType = Image\n"); + fprintf(FMHD, "NDims = 3\n"); + fprintf(FMHD, "DimSize = %d %d %d\n", dimx, dimy, dimz); + fprintf(FMHD, "ElementType = MET_SHORT\n"); + fprintf(FMHD, "ElementSpacing = %.1f %.1f %.1f\n", resolution, + resolution, resolution); + fprintf(FMHD, "ElementByteOrderMSB = False\n"); + fprintf(FMHD, "ElementDataFile = %s\n", "invasion_diameters.raw"); + fprintf(FMHD, "HeaderSize = 0\n"); + fclose(FMHD); + } + + FILE *log_file = fopen("injection_output.csv", "a"); + fprintf(log_file, "%d %d %f %d %f %d %f\n", step, D, D * resolution, + injectedVolume, injectedVolume / (1.0 * NP), displacedVolume, + displacedVolume / (1.0 * NP)); + fclose(log_file); + + return D; +} + +int main(int argc, char *argv[]) { + + if (argc != 2) + ERROR("Wrong number of parameters."); + + Full_Morphology fm(argc, argv); + + const int nsteps = fm.diameter.size(); + for (int step = 0; step < nsteps; step++) { + + int d = fm.calc(step); + cout << "Step " << step << ", D = " << d << " px." << endl; + } + + return 0; +} diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index cafa7a52..977d304b 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -66,7 +66,8 @@ int main( int argc, char **argv ) // structure and allocate variables ColorModel.Initialize(); // initializing the model will set initial conditions for variables - if (SimulationMode == "legacy"){ + auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "default" ); + if (PROTOCOL == "sw_steady") { ColorModel.Run(); } else { @@ -75,7 +76,6 @@ int main( int argc, char **argv ) bool ContinueSimulation = true; /* Variables for simulation protocols */ - auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "default" ); /* image sequence protocol */ int IMAGE_INDEX = 0; int IMAGE_COUNT = 0; diff --git a/tests/test_MPI.cpp b/tests/test_MPI.cpp index c6d12011..f14ba219 100644 --- a/tests/test_MPI.cpp +++ b/tests/test_MPI.cpp @@ -1652,7 +1652,9 @@ int main( int argc, char *argv[] ) } // Test GPU aware MPI +#ifdef USE_CUDA test_GPU_aware( &ut ); +#endif } // Limit the scope so objects are destroyed