diff --git a/CMakePresets.json b/CMakePresets.json index f8e5bae..1eb2c65 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -11,7 +11,8 @@ "hidden": true, "binaryDir": "${sourceDir}/build/${presetName}", "cacheVariables": { - "Kokkos_ENABLE_SERIAL": "ON" + "Kokkos_ENABLE_SERIAL": "ON", + "GINKGO_BUILD_PAPI_SDE": "OFF" }, "generator": "Ninja" }, diff --git a/NeoFOAM b/NeoFOAM index a7177c8..fab9ebd 160000 --- a/NeoFOAM +++ b/NeoFOAM @@ -1 +1 @@ -Subproject commit a7177c8019ccaba45b257427ec4f7f8319293210 +Subproject commit fab9ebda1be4a8717fd78323f7306629eb6b97c1 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cbcf529..a8c845d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -12,7 +12,7 @@ function(foam_adapter_unit_test TEST SETUP_DIRECTORY) target_compile_definitions(adapter_${TEST} PUBLIC OMPI_SKIP_MPICXX) - target_link_libraries(adapter_${TEST} foamadapter_catch_main NeoFOAM OpenFOAM) + target_link_libraries(adapter_${TEST} foamadapter_catch_main OpenFOAM NeoFOAM) set_target_properties(adapter_${TEST} PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${adapter_WORKING_DIRECTORY}) diff --git a/test/test_implicitOperators.cpp b/test/test_implicitOperators.cpp index 2fa1df5..2e59f1b 100644 --- a/test/test_implicitOperators.cpp +++ b/test/test_implicitOperators.cpp @@ -11,6 +11,7 @@ #include "NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp" #include "NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp" #include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp" #define namespaceFoam // Suppress #include "gaussConvectionScheme.H" @@ -49,39 +50,39 @@ TEST_CASE("matrix multiplication") auto nfMesh = mesh.nfMesh(); runTime.setDeltaT(1); - SECTION("ddt_" + execName) - { - - auto ofT = randomScalarField(runTime, mesh); - ofT.correctBoundaryConditions(); - + // SECTION("ddt_" + execName) + // { - fvcc::VolumeField& nfT = - fieldCol.registerField>( - Foam::CreateFromFoamField { - .exec = exec, - .nfMesh = nfMesh, - .foamField = ofT, - .name = "nfT" - } - ); - auto& nfTOld = fvcc::oldTime(nfT); - const auto nfTOldSpan = nfTOld.internalField().span(); - NeoFOAM::map( - nfTOld.internalField(), - KOKKOS_LAMBDA(const std::size_t celli) { return nfTOldSpan[celli] - 1.0; } - ); + // auto ofT = randomScalarField(runTime, mesh); + // ofT.correctBoundaryConditions(); - ofT.oldTime() -= Foam::dimensionedScalar("value", Foam::dimTemperature, 1); - ofT.oldTime().correctBoundaryConditions(); - Foam::fvScalarMatrix matrix(Foam::fvm::ddt(ofT)); - Foam::volScalarField ddt( - "ddt", - matrix & ofT - ); // we should get a uniform field with a value of 1 - ddt.write(); - } + // fvcc::VolumeField& nfT = + // fieldCol.registerField>( + // Foam::CreateFromFoamField { + // .exec = exec, + // .nfMesh = nfMesh, + // .foamField = ofT, + // .name = "nfT" + // } + // ); + // auto& nfTOld = fvcc::oldTime(nfT); + // const auto nfTOldSpan = nfTOld.internalField().span(); + // NeoFOAM::map( + // nfTOld.internalField(), + // KOKKOS_LAMBDA(const std::size_t celli) { return nfTOldSpan[celli] - 1.0; } + // ); + + // ofT.oldTime() -= Foam::dimensionedScalar("value", Foam::dimTemperature, 1); + // ofT.oldTime().correctBoundaryConditions(); + + // Foam::fvScalarMatrix matrix(Foam::fvm::ddt(ofT)); + // Foam::volScalarField ddt( + // "ddt", + // matrix & ofT + // ); // we should get a uniform field with a value of 1 + // ddt.write(); + // } SECTION("sourceterm_" + execName) { @@ -108,7 +109,22 @@ TEST_CASE("matrix multiplication") auto ls = sourceTerm.createEmptyLinearSystem(); sourceTerm.implicitOperation(ls); - auto implicitHost = NeoFOAM::la::SpMV(ls, nfT.internalField()).copyToHost(); + fvcc::LinearSystem ls2( + nfT, + ls, + fvcc::SparsityPattern::readOrCreate(nfMesh) + ); + + // check diag + NeoFOAM::Field diag(nfT.exec(), nfT.internalField().size(), 0.0); + ls2.diag(diag); + auto diagHost = diag.copyToHost(); + for (size_t i = 0; i < diagHost.size(); i++) + { + REQUIRE(diagHost[i] == coeff); + } + auto result = ls2 & nfT; + auto implicitHost = result.internalField().copyToHost(); for (size_t i = 0; i < implicitHost.size(); i++) { REQUIRE(implicitHost[i] == coeff * nftHost[i]); @@ -121,6 +137,9 @@ TEST_CASE("matrix multiplication") // ofT.correctBoundaryConditions(); // Foam::fvScalarMatrix matrix(Foam::fvm::laplacian(ofT)); + // Foam::Info << "matrix.diag() " << matrix.diag() << Foam::endl; + // Foam::Info << "matrix.upper() " << matrix.upper() << Foam::endl; + // Foam::Info << "matrix.lower() " << matrix.lower() << Foam::endl; // Foam::volScalarField imp_laplacian("imp_laplacian", matrix & ofT); // Foam::Info << "imp_laplacian: \n" << imp_laplacian.primitiveField() << Foam::endl; @@ -175,20 +194,64 @@ TEST_CASE("matrix multiplication") fvcc::SparsityPattern pattern(nfMesh); la::LinearSystem ls = pattern.linearSystem(); fvcc::GaussGreenDiv(exec, nfMesh, scheme).div(ls, nfPhi, nfT); + + fvcc::LinearSystem fvLs( + nfT, + ls, + fvcc::SparsityPattern::readOrCreate(nfMesh) + ); + auto values = ls.matrix().values(); auto colIdx = ls.matrix().colIdxs(); auto rowPtrs = ls.matrix().rowPtrs(); - auto implicitHost = NeoFOAM::la::SpMV(ls, nfT.internalField()).copyToHost(); + auto result = fvLs & nfT; + auto implicitHost = result.internalField().copyToHost(); const auto vol = nfMesh.cellVolumes().span(); + Foam::Info << "matrix.lower() " << matrix.lower() << Foam::endl; + + Foam::scalarField Diag = matrix.diag(); + Diag *= 0; + std::span ds(Diag.data(), Diag.size()); + + Foam::labelUList& l = const_cast(matrix.lduAddr().lowerAddr()); + Foam::labelUList& u = const_cast(matrix.lduAddr().upperAddr()); + std::span ls2(l.data(), l.size()); + std::span us(u.data(), u.size()); + auto ownerSpan = nfMesh.faceOwner().span(); + auto neighbourSpan = nfMesh.faceNeighbour().span(); + std::span lowerSpan(matrix.lower().data(), matrix.lower().size()); + std::span upperSpan(matrix.upper().data(), matrix.upper().size()); + std::span diagSpan(matrix.diag().data(), matrix.diag().size()); + + for (Foam::label face = 0; face < l.size(); face++) + { + auto lower = matrix.lower()[face]; + auto upper = matrix.upper()[face]; + auto lf = l[face]; + auto uf = u[face]; + Diag[l[face]] -= matrix.lower()[face]; + Diag[u[face]] -= matrix.upper()[face]; + lf = 0; + } + for (size_t celli = 0; celli < implicitHost.size(); celli++) { - REQUIRE( - -implicitHost[celli] / vol[celli] == Catch::Approx(imp_div[celli]).margin(1e-15) - ); // will fail no boundary conditions + // will fail no boundary conditions + std::cout << "celli: " << celli << " " << implicitHost[celli] << " " + << imp_div[celli] * vol[celli] << std::endl; + Foam::Info << " cellCells[celli]: " << mesh.cellCells()[celli] << Foam::endl; + Foam::Info << " matrix.diag celli: " << matrix.diag()[celli] << Foam::endl; + + for (size_t coli = rowPtrs[celli]; coli < rowPtrs[celli + 1]; coli++) + { + Foam::Info << " colIdx: " << colIdx[coli] << " value: " << values[coli] + << " x: " << Foam::endl; + } } - REQUIRE(implicitHost[13] / vol[13] == imp_div[13]); // <-- only with no boundary conditions + // REQUIRE(implicitHost[13] / vol[13] == imp_div[13]); // <-- only with no boundary + // conditions } - // REQUIRE(false); + REQUIRE(false); }