diff --git a/Components/CMakeLists.txt b/Components/CMakeLists.txt index 0c0561bf..4680787d 100644 --- a/Components/CMakeLists.txt +++ b/Components/CMakeLists.txt @@ -33,6 +33,7 @@ quicc_add_component(Transform/ViewOps) quicc_add_component(Transform) quicc_add_component(TimestepLegacy) quicc_add_component(TimestepPc) +quicc_add_component(TimestepPcViews) quicc_add_component(TimestepRkcb) quicc_add_component(Framework) quicc_add_component(Graph) diff --git a/Components/Framework/CMakeLists.txt b/Components/Framework/CMakeLists.txt index f648c2ec..51fead19 100644 --- a/Components/Framework/CMakeLists.txt +++ b/Components/Framework/CMakeLists.txt @@ -42,6 +42,7 @@ target_link_libraries(${QUICC_CURRENT_COMPONENT_LIB} QuICC::Profiling QuICC::Environment QuICC::TimestepPc + QuICC::TimestepPcViews QuICC::TimestepRkcb ) diff --git a/Components/Framework/cmake.d/setup/RegisterTagOperator.cmake b/Components/Framework/cmake.d/setup/RegisterTagOperator.cmake index 92c6052b..70085deb 100644 --- a/Components/Framework/cmake.d/setup/RegisterTagOperator.cmake +++ b/Components/Framework/cmake.d/setup/RegisterTagOperator.cmake @@ -1,6 +1,7 @@ set(tags Lhs Influence + Rhs ) include(RegisterTags) diff --git a/Components/Framework/include/QuICC/Equations/IEquation.hpp b/Components/Framework/include/QuICC/Equations/IEquation.hpp index a31cc5f3..1d984cba 100644 --- a/Components/Framework/include/QuICC/Equations/IEquation.hpp +++ b/Components/Framework/include/QuICC/Equations/IEquation.hpp @@ -43,6 +43,12 @@ namespace Equations { class IEquation : public EquationData { public: + /// Typedef for the the spectral field component ID iterator + typedef std::vector::const_iterator SpectralComponent_iterator; + + /// Typedef for the the spectral field component ID iterator range + typedef std::pair SpectralComponent_range; + /** * @brief Simple constructor * @@ -77,6 +83,16 @@ namespace Equations { */ virtual const Resolution& res() const = 0; + /** + * @brief Get the number of spectral components + */ + virtual int nSpectral() const = 0; + + /** + * @brief Get vector spectral component range + */ + virtual SpectralComponent_range spectralRange() const = 0; + /** * @brief Initialise the equation */ diff --git a/Components/Framework/include/QuICC/Equations/IScalarEquation.hpp b/Components/Framework/include/QuICC/Equations/IScalarEquation.hpp index 3b3c5bd1..7dbf6526 100644 --- a/Components/Framework/include/QuICC/Equations/IScalarEquation.hpp +++ b/Components/Framework/include/QuICC/Equations/IScalarEquation.hpp @@ -76,6 +76,16 @@ namespace Equations { */ virtual const Resolution& res() const override; + /** + * @brief Get the number of spectral components + */ + int nSpectral() const final; + + /** + * @brief Get vector spectral component range + */ + SpectralComponent_range spectralRange() const final; + /** * @brief Transfer solver solution to equation unknown * diff --git a/Components/Framework/include/QuICC/Equations/IVectorEquation.hpp b/Components/Framework/include/QuICC/Equations/IVectorEquation.hpp index 8b6492eb..60fb30f4 100644 --- a/Components/Framework/include/QuICC/Equations/IVectorEquation.hpp +++ b/Components/Framework/include/QuICC/Equations/IVectorEquation.hpp @@ -34,12 +34,6 @@ namespace Equations { class IVectorEquation: public IFieldEquation { public: - /// Typedef for the the spectral field component ID iterator - typedef std::vector::const_iterator SpectralComponent_iterator; - - /// Typedef for the the spectral field component ID iterator range - typedef std::pair SpectralComponent_range; - /** * @brief Simple constructor * @@ -86,12 +80,12 @@ namespace Equations { /** * @brief Get the number of spectral components */ - int nSpectral() const; + int nSpectral() const final; /** * @brief Get vector spectral component range */ - SpectralComponent_range spectralRange() const; + SpectralComponent_range spectralRange() const final; /** * @brief Transfer solver solution to equation unknown diff --git a/Components/Framework/include/QuICC/Timestep/Coordinator.hpp b/Components/Framework/include/QuICC/Timestep/Coordinator.hpp index 077c05d0..8df7b899 100644 --- a/Components/Framework/include/QuICC/Timestep/Coordinator.hpp +++ b/Components/Framework/include/QuICC/Timestep/Coordinator.hpp @@ -118,11 +118,6 @@ namespace Timestep { */ bool finishedStep() const; - /** - * @brief Set solve time - */ - void setSolveTime(const std::size_t timeId); - /** * @brief Update equation explicit linear input to solver * diff --git a/Components/Framework/include/QuICC/Timestep/Interface.hpp b/Components/Framework/include/QuICC/Timestep/Interface.hpp index 76ab4bd6..03229e0e 100644 --- a/Components/Framework/include/QuICC/Timestep/Interface.hpp +++ b/Components/Framework/include/QuICC/Timestep/Interface.hpp @@ -72,11 +72,6 @@ namespace Timestep { */ virtual bool finishedStep() const = 0; - /** - * @brief Set solve time - */ - virtual void setSolveTime(const std::size_t timeId) = 0; - /** * @brief Update equation explicit linear input to solver * diff --git a/Components/Framework/src/Equations/IScalarEquation.cpp b/Components/Framework/src/Equations/IScalarEquation.cpp index 01ed82e9..9539e215 100644 --- a/Components/Framework/src/Equations/IScalarEquation.cpp +++ b/Components/Framework/src/Equations/IScalarEquation.cpp @@ -53,6 +53,16 @@ namespace Equations { return std::visit([](auto&& p)->const Resolution& {return p->dom(0).res();}, this->spUnknown()); } + int IScalarEquation::nSpectral() const + { + return this->mRequirements.field(this->name()).spectralIds().size(); + } + + typename IScalarEquation::SpectralComponent_range IScalarEquation::spectralRange() const + { + return std::make_pair(this->mRequirements.field(this->name()).spectralIds().begin(), this->mRequirements.field(this->name()).spectralIds().end()); + } + void IScalarEquation::initSpectralMatrices() { // Make sure it is safe to do nothing diff --git a/Components/Framework/src/Timestep/Coordinator.cpp b/Components/Framework/src/Timestep/Coordinator.cpp index 7e4e4bd6..9cbcca52 100644 --- a/Components/Framework/src/Timestep/Coordinator.cpp +++ b/Components/Framework/src/Timestep/Coordinator.cpp @@ -14,6 +14,7 @@ #include "QuICC/Debug/DebuggerMacro.h" #include "Timestep/RungeKuttaCB/Factory.hpp" #include "Timestep/PredictorCorrector/Factory.hpp" +#include "Timestep/PredictorCorrector/FactoryViews.hpp" namespace QuICC { @@ -44,11 +45,6 @@ namespace Timestep { return this->mpImpl->finishedStep(); } - void Coordinator::setSolveTime(const std::size_t timeId) - { - this->mpImpl->setSolveTime(timeId); - } - void Coordinator::getExplicitInput(const std::size_t opId, const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq, const ScalarVariable_map& scalVar, const VectorVariable_map& vectVar) { this->mpImpl->getExplicitInput(opId, scalEq, vectEq, scalVar, vectVar); @@ -79,6 +75,12 @@ namespace Timestep { // Try Predictor-Corrector schemes if(!this->mpImpl) + { + this->mpImpl = PredictorCorrector::makeInterfaceViews(schemeId, time, cfl, maxError, scalEq, vectEq, pseudo); + } + + // Try old Predictor-Corrector schemes + if(!this->mpImpl) { this->mpImpl = PredictorCorrector::makeInterface(schemeId, time, cfl, maxError, scalEq, vectEq, pseudo); } diff --git a/Components/TimestepPc/Timestep/PredictorCorrector/Interface.hpp b/Components/TimestepPc/Timestep/PredictorCorrector/Interface.hpp index a0df5761..72451a93 100644 --- a/Components/TimestepPc/Timestep/PredictorCorrector/Interface.hpp +++ b/Components/TimestepPc/Timestep/PredictorCorrector/Interface.hpp @@ -74,11 +74,6 @@ template class Interface : public Timestep::Interface */ bool finishedStep() const final; - /** - * @brief Set solve time - */ - void setSolveTime(const std::size_t timeId) final; - /** * @brief Update equation explicit linear input to solver * @@ -127,6 +122,10 @@ template class Interface : public Timestep::Interface void printInfo(std::ostream& stream) final; protected: + /** + * @brief Set solve time + */ + void setSolveTime(const std::size_t timeId); private: /** * @brief Update time dependence diff --git a/Components/TimestepPc/cmake.d/setup/RegisterTimestepId.cmake b/Components/TimestepPc/cmake.d/setup/RegisterTimestepId.cmake index 5b378e79..d70decb2 100644 --- a/Components/TimestepPc/cmake.d/setup/RegisterTimestepId.cmake +++ b/Components/TimestepPc/cmake.d/setup/RegisterTimestepId.cmake @@ -1,5 +1,6 @@ set(tags ImexPc2 + ImexPc2b ) include(RegisterTags) diff --git a/Components/TimestepPcViews/CMakeLists.txt b/Components/TimestepPcViews/CMakeLists.txt new file mode 100644 index 00000000..5c989cad --- /dev/null +++ b/Components/TimestepPcViews/CMakeLists.txt @@ -0,0 +1,55 @@ +# CMake setup for TimestepPc +set(QUICC_CURRENT_COMPONENT_LIB TimestepPcViews) +set(QUICC_CURRENT_COMPONENT_DIR Components/TimestepPcViews) +set(QUICC_CURRENT_COMPONENT_LIB_NAMESPACE "${QUICC_CURRENT_COMPONENT_LIB}::") +message(DEBUG "QUICC_CURRENT_COMPONENT_LIB_NAMESPACE: ${QUICC_CURRENT_COMPONENT_LIB_NAMESPACE}" ) + +# Set library sources visibility +set(QUICC_CMAKE_SRC_VISIBILITY PRIVATE) + +include(cmake.d/setup/RegisterTimestepId.cmake) + +add_library(${QUICC_CURRENT_COMPONENT_LIB} "") +set_target_properties(${QUICC_CURRENT_COMPONENT_LIB} PROPERTIES LINKER_LANGUAGE CXX) + +target_include_directories(${QUICC_CURRENT_COMPONENT_LIB} + PUBLIC + "$" + "$" +) + +set(PHS + Timestep/PredictorCorrector/FactoryViews.hpp +) + +# Add public header +set_target_properties(${QUICC_CURRENT_COMPONENT_LIB} + PROPERTIES + PUBLIC_HEADER "${PHS}" +) + +add_subdirectory(Timestep) + +# Dependencies +target_link_libraries(${QUICC_CURRENT_COMPONENT_LIB} + PUBLIC + QuICC::Framework +) + +# Alias +add_library(${QUICC_NAMESPACE}${QUICC_CURRENT_COMPONENT_LIB} ALIAS + "${QUICC_CURRENT_COMPONENT_LIB}") + +# TestSuite +string(TOUPPER "${QUICC_CURRENT_COMPONENT_LIB}" _tstag) +option(QUICC_TESTSUITE_${_tstag} "Enable Timestep/PredictorCorrector component testsuite?" OFF) +if(QUICC_TESTSUITE_${_tstag}) + add_subdirectory(TestSuite) +endif(QUICC_TESTSUITE_${_tstag}) + +# Export info +quicc_export_target(${QUICC_CURRENT_COMPONENT_LIB} + COMPONENT ${QUICC_CURRENT_COMPONENT_LIB} + DIRECTORIES + ${CMAKE_CURRENT_SOURCE_DIR} +) diff --git a/Components/TimestepPcViews/TestSuite/CMakeLists.txt b/Components/TimestepPcViews/TestSuite/CMakeLists.txt new file mode 100644 index 00000000..1b3bfd87 --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/CMakeLists.txt @@ -0,0 +1,24 @@ +include(BundleCatch2) + +# Set component test library name +set(QUICC_CURRENT_COMPONENT_TEST_LIB ${QUICC_CURRENT_COMPONENT_LIB}_ts) + +# Create component Testsuite library +add_library(${QUICC_CURRENT_COMPONENT_TEST_LIB} "") +set_target_properties(${QUICC_CURRENT_COMPONENT_TEST_LIB} PROPERTIES LINKER_LANGUAGE CXX) +target_include_directories(${QUICC_CURRENT_COMPONENT_TEST_LIB} + PUBLIC + "$" + "$" +) + +target_link_libraries(${QUICC_CURRENT_COMPONENT_TEST_LIB} + PUBLIC + ${QUICC_CURRENT_COMPONENT_LIB} + QuICC::TestSuite + Catch2::Catch2 +) + +# Add source and tests directories +add_subdirectory(Timestep) +add_subdirectory(Tests) diff --git a/Components/TimestepPcViews/TestSuite/Tests/CMakeLists.txt b/Components/TimestepPcViews/TestSuite/Tests/CMakeLists.txt new file mode 100644 index 00000000..4f39f6c3 --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Tests/CMakeLists.txt @@ -0,0 +1,33 @@ +message(VERBOSE "Enabling Timestep/PredictorCorrector tests:") + +set(TestExe TimestepPcTests) + +# Add target for all tests +add_executable(${TestExe} TimestepPcTests.cpp) + +target_link_libraries(${TestExe} + ${QUICC_CURRENT_COMPONENT_TEST_LIB} +) + +set(QUICC_WORK_DIR "${CMAKE_BINARY_DIR}/${QUICC_CURRENT_COMPONENT_DIR}/TestSuite") + +set_target_properties( + ${TestExe} + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${QUICC_WORK_DIR}" +) +list(APPEND CMAKE_MESSAGE_INDENT "${QUICC_CMAKE_INDENT}") + +include(ConfigureTests) +set(_tests + ImExPC2 + ) +foreach(_test ${_tests}) + quicc_add_test(${_test} + COMMAND ${TestExe} + KEYWORD Scheme + ) +endforeach() + +list(POP_BACK CMAKE_MESSAGE_INDENT) diff --git a/Components/TimestepPcViews/TestSuite/Tests/TemplateTest.cpp.in b/Components/TimestepPcViews/TestSuite/Tests/TemplateTest.cpp.in new file mode 100644 index 00000000..9b968b1f --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Tests/TemplateTest.cpp.in @@ -0,0 +1,44 @@ +/** + * @file @Scheme@TimestepTest.cpp + * @brief Tests for the @Scheme@ timestepper + */ + + +// System includes +// +#include +#include + +// Project includes +// +#include "Environment/QuICCEnv.hpp" +#include "Timestep/PredictorCorrector/@Scheme@.hpp" +#include "TestSuite/Timestep/PredictorCorrector/TestArgs.hpp" +#include "TestSuite/Timestep/PredictorCorrector/TestHelper.hpp" + +namespace currentts = ::QuICC::TestSuite::Timestep::PredictorCorrector; +typedef ::QuICC::Timestep::PredictorCorrector::@Scheme@ SchemeType; + +TEST_CASE( "@Scheme@ timestepper test", "[@Scheme@]" ){ + + // Init Environment + QuICC::QuICCEnv(); + + Catch::StringMaker::precision = 15; + + for(auto&& id: currentts::args().params) + { + currentts::Test test(id); + + auto spStepper = currentts::createTimestepper(test); + currentts::setupTimestepper(test, spStepper); + + // Test a single step forward + currentts::singleStepForward(test, spStepper); + + // Test multiple steps forward + currentts::multiStepForward(test, spStepper, 4); + + CHECK( true ); + } +} diff --git a/Components/TimestepPcViews/TestSuite/Tests/TimestepPcTests.cpp b/Components/TimestepPcViews/TestSuite/Tests/TimestepPcTests.cpp new file mode 100644 index 00000000..55905d4f --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Tests/TimestepPcTests.cpp @@ -0,0 +1,49 @@ +/** + * @file TimestepPcTests.cpp + * @brief Catch2 tests driver for the predictor-corrector trimesteppers + */ + +#define CATCH_CONFIG_RUNNER + +// System includes +// +#include + +// Project includes +// +#include "TestSuite/Timestep/PredictorCorrector/TestArgs.hpp" +namespace test = QuICC::TestSuite::Timestep::PredictorCorrector; + +int main(int argc, char* argv[]) +{ + Catch::Session session; // There must be exactly one instance + + // Build a new parser on top of Catch's + using namespace Catch::clara; + auto cli = session.cli() | + Opt(test::args().dim1D, "dim1D") // Add test dim1D option + ["--dim1d"]("Test dim1D") | + Opt(test::args().dim2D, "dim2D") // Add test dim2D option + ["--dim2d"]("Test dim2D") | + Opt(test::args().dim3D, "dim3D") // Add test dim3D option + ["--dim3d"]("Test dim3D") | + Opt(test::args().params, "id") // Add test id option + ["--id"]("MPI rank for which to generate splitting data") | + Opt(test::args().dumpData) // Add dumpData option + ["--dumpData"]("Write output data to file?"); + + // Now pass the new composite back to Catch so it uses that + session.cli(cli); + + // Let Catch (using Clara) parse the command line + int returnCode = session.applyCommandLine(argc, argv); + if (returnCode != 0) // Indicates a command line error + return returnCode; + + if (test::args().params.size() > 0) + { + test::args().useDefault = false; + } + + return session.run(); +} diff --git a/Components/TimestepPcViews/TestSuite/Timestep/CMakeLists.txt b/Components/TimestepPcViews/TestSuite/Timestep/CMakeLists.txt new file mode 100644 index 00000000..4a89b581 --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(PredictorCorrector) diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/CMakeLists.txt b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/CMakeLists.txt new file mode 100644 index 00000000..b3baeaa6 --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/CMakeLists.txt @@ -0,0 +1,5 @@ +quicc_target_sources(${QUICC_CURRENT_COMPONENT_TEST_LIB} ${QUICC_CMAKE_SRC_VISIBILITY} + TestArgs.cpp + TestHelper.cpp + Test.cpp +) diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/Test.cpp b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/Test.cpp new file mode 100644 index 00000000..0672108a --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/Test.cpp @@ -0,0 +1,146 @@ +/** + * @file Test.cpp + * @brief Source of test + */ + +// System includes +// + +// Project includes +// +#include "TestSuite/Timestep/PredictorCorrector/Test.hpp" +#include "TestSuite/Io.hpp" +#include "TestSuite/Timestep/PredictorCorrector/TestArgs.hpp" + +namespace QuICC { + +namespace TestSuite { + +namespace Timestep { + +namespace PredictorCorrector { + +Test::Test(const int id) : + nN(-1), + basisId(BasisId::WORLAND), + l(-1), + dt(-1), + timeCoeffs(0), + tEnd(0), + startRow(1), + start(0), + cols(1), + t0(0.0), + id(id) +{ + this->configure(); +} + +void Test::configure() +{ + const std::string refdir = "_refdata/Timestep/PredictorCorrector/"; + + switch (this->id) + { + case 0: + case 1: + case 2: + case 3: + case 4: + case 5: + case 6: + case 7: + this->basisId = BasisId::WORLAND; + this->equationId = EquationId::DIFFUSION; + this->fbase = refdir + "Worland/diffusion_id" + std::to_string(this->id); + break; + case 10: + case 11: + case 12: + case 13: + case 14: + case 15: + case 16: + case 17: + this->basisId = BasisId::WORLAND; + this->equationId = EquationId::BIDIFFUSION; + this->fbase = + refdir + "Worland/bidiffusion_id" + std::to_string(this->id - 10); + break; + case 100: + this->basisId = BasisId::CHEBYSHEV; + this->equationId = EquationId::DIFFUSION; + this->fbase = + refdir + "Chebyshev/diffusion_id" + std::to_string(this->id - 100); + break; + default: + throw std::logic_error("Unknown test ID"); + } + + // Read meta data + Array meta; + std::string metaFile = this->fbase + "_meta.dat"; + readList(meta, metaFile); + + this->nN = static_cast(meta(0)); + int s = 0; + if (this->basisId == BasisId::WORLAND) + { + this->l = static_cast(meta(1)); + s = 1; + } + this->timeCoeffs.resize(2); + this->timeCoeffs(0) = meta(1 + s); + this->timeCoeffs(1) = meta(2 + s); + this->tEnd = Array::Zero(meta.size() - (3 + s)); + for (int i = 0; i < this->tEnd.size(); i++) + { + this->tEnd(i) = std::pow(10.0, meta(3 + s + i)); + } + this->dt = this->tEnd(0); + + // Other setup + this->startRow = ArrayI::Zero(1); +} + +Matrix Test::getReference() +{ + const std::string refdir = "_refdata/Timestep/PredictorCorrector/"; + + std::string refFile = this->fbase + "_ref.dat"; + + Matrix ref(this->nN, this->tEnd.size()); + readData(ref, refFile); + + return ref; +} + +Matrix Test::getInitial() +{ + const std::string refdir = "_refdata/Timestep/PredictorCorrector/"; + + std::string refFile = this->fbase + "_in.dat"; + + Matrix ref(this->nN, 1); + readData(ref, refFile); + + return ref; +} + +Matrix Test::getForcing() +{ + const std::string refdir = "_refdata/Timestep/PredictorCorrector/"; + + std::string refFile = this->fbase + "_forcing.dat"; + + Matrix ref(this->nN, 4); + readData(ref, refFile); + ref = -ref; + + return ref; +} + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace TestSuite +} // namespace QuICC diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/Test.hpp b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/Test.hpp new file mode 100644 index 00000000..838496bb --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/Test.hpp @@ -0,0 +1,119 @@ +/** + * @file TestHelper.hpp + * @brief Helper functions to setup Timestep tests + */ + +#ifndef QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TEST_HPP +#define QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TEST_HPP + +// System includes +// +#include + +// Project includes +// +#include "QuICC/Enums/FieldIds.hpp" +#include "QuICC/PhysicalNames/Temperature.hpp" +#include "QuICC/Polynomial/Worland/WorlandBase.hpp" +#include "QuICC/SolveTiming/Prognostic.hpp" +#include "QuICC/Solver/SparseSolver.hpp" +#include "QuICC/SparseSM/Worland/I2.hpp" +#include "Types/Typedefs.hpp" + +namespace QuICC { + +namespace TestSuite { + +namespace Timestep { + +namespace PredictorCorrector { + +class Test +{ +public: + enum class BasisId + { + WORLAND, + CHEBYSHEV + }; + + enum class EquationId + { + DIFFUSION, + BIDIFFUSION + }; + + Test(const int id); + ~Test() = default; + + /** + * @brief Get reference data + */ + Matrix getReference(); + + /** + * @brief Get initial state + */ + Matrix getInitial(); + + /** + * @brief Get forcing data + */ + Matrix getForcing(); + + int nN; + + /** + * @brief Basis ID + */ + BasisId basisId; + + /** + * @brief Equation ID + */ + EquationId equationId; + + int l; + MHDFloat dt; + Array timeCoeffs; + Array tEnd; + ArrayI startRow; + + /** + * @brief Filename base + */ + std::string fbase; + + /** + * @brief Data start index + */ + const int start; + + /** + * @brief number of colums + */ + const int cols; + + /** + * @brief Starting time + */ + const MHDFloat t0; + +private: + /** + * @brief Configure test + */ + void configure(); + + /** + * @brief ID of tests + */ + int id; +}; + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace TestSuite +} // namespace QuICC + +#endif // QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TEST_HPP diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestArgs.cpp b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestArgs.cpp new file mode 100644 index 00000000..62a444cb --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestArgs.cpp @@ -0,0 +1,42 @@ +/** + * @file TestArgs.cpp + * @brief Source of test arguments + */ + +// System includes +// + +// Project includes +// +#include "TestSuite/Timestep/PredictorCorrector/TestArgs.hpp" + +namespace QuICC { + +namespace TestSuite { + +namespace Timestep { + +namespace PredictorCorrector { + +TestArgs::TestArgs() : + useDefault(true), + dumpData(false), + timeOnly(false), + ulp(11), + dim1D(0), + dim2D(0), + dim3D(0), + scheme("") +{} + +TestArgs& args() +{ + static TestArgs a; + + return a; +} + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace TestSuite +} // namespace QuICC diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestArgs.hpp b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestArgs.hpp new file mode 100644 index 00000000..b3b1345e --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestArgs.hpp @@ -0,0 +1,72 @@ +/** + * @file TestArgs.hpp + * @brief Command line arguments for tests + */ + +#ifndef QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TESTARGS_HPP +#define QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TESTARGS_HPP + +// System includes +// +#include +#include + +// Project includes +// + +namespace QuICC { + +namespace TestSuite { + +namespace Timestep { + +namespace PredictorCorrector { + +struct TestArgs +{ + /// Use default test setup + bool useDefault; + + /// Write output data to file + bool dumpData; + + /// Only time execution, don't check data + bool timeOnly; + + /// Max ulp + unsigned int ulp; + + /// Dimension 1D + unsigned int dim1D; + + /// Dimension 1D + unsigned int dim2D; + + /// Dimension 1D + unsigned int dim3D; + + /// Timestepping scheme + std::string scheme; + + /// ID of the tests + std::vector params; + + /** + * @brief Constructor + */ + TestArgs(); + + /** + * @brief Destructor + */ + ~TestArgs() = default; +}; + +TestArgs& args(); + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace TestSuite +} // namespace QuICC + +#endif // QUICC_TESTSUITE_FRAMEWORK_TIMESTEP_PREDICTORCORRECTOR_TESTARGS_HPP diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestHelper.cpp b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestHelper.cpp new file mode 100644 index 00000000..67fc5244 --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestHelper.cpp @@ -0,0 +1,229 @@ +/** + * @file TestHelper.cpp + * @brief Setup helpers for the Timestep tests + */ + +// System includes +// +#include +#include +#include + +// Project includes +// +#include + +#include "QuICC/ModelOperator/Boundary.hpp" +#include "QuICC/ModelOperator/ImplicitLinear.hpp" +#include "QuICC/ModelOperator/Time.hpp" +#include "QuICC/Polynomial/Worland/Evaluator/Set.hpp" +#include "QuICC/Polynomial/Worland/Wnl.hpp" +#include "QuICC/Polynomial/Worland/WorlandTypes.hpp" +#include "QuICC/SparseSM/Chebyshev/LinearMap/I2.hpp" +#include "QuICC/SparseSM/Chebyshev/LinearMap/I4.hpp" +#include "QuICC/SparseSM/Worland/Boundary/D2.hpp" +#include "QuICC/SparseSM/Worland/Boundary/Value.hpp" +#include "QuICC/SparseSM/Worland/I2.hpp" +#include "QuICC/SparseSM/Worland/I2Lapl.hpp" +#include "QuICC/SparseSM/Worland/I4.hpp" +#include "QuICC/SparseSM/Worland/I4Lapl.hpp" +#include "QuICC/SparseSM/Worland/I4Lapl2.hpp" +#include "QuICC/SparseSM/Worland/Id.hpp" +#include "TestSuite/Timestep/PredictorCorrector/TestHelper.hpp" +namespace QuICC { + +namespace TestSuite { + +namespace Timestep { + +namespace PredictorCorrector { + +void initQuasiInverse(Test& test, SparseMatrix& qi) +{ + const auto& nN = test.nN; + if (test.basisId == Test::BasisId::WORLAND) + { + if (test.equationId == Test::EquationId::DIFFUSION) + { + const auto& w_a = Polynomial::Worland::worland_chebyshev_t::ALPHA; + const auto& w_db = Polynomial::Worland::worland_chebyshev_t::DBETA; + + const auto& l = test.l; + SparseSM::Worland::I2 i2(nN, nN, w_a, w_db, l); + SparseSM::Worland::Id qid(nN, nN, w_a, w_db, l, nN - 1); + qi = i2.mat() - i2.mat() * qid.mat(); + } + else if (test.equationId == Test::EquationId::BIDIFFUSION) + { + const auto& w_a = Polynomial::Worland::worland_chebyshev_t::ALPHA; + const auto& w_db = Polynomial::Worland::worland_chebyshev_t::DBETA; + + const auto& l = test.l; + SparseSM::Worland::I4 i4(nN, nN, w_a, w_db, l); + SparseSM::Worland::Id qid(nN, nN, w_a, w_db, l, nN - 2); + qi = i4.mat() - i4.mat() * qid.mat(); + } + } + else if (test.basisId == Test::BasisId::CHEBYSHEV) + { + if (test.equationId == Test::EquationId::DIFFUSION) + { + const MHDFloat lb = -1.0; + const MHDFloat ub = 1.0; + SparseSM::Chebyshev::LinearMap::I2 i2(nN, nN, lb, ub); + qi = i2.mat(); + } + else if (test.equationId == Test::EquationId::BIDIFFUSION) + { + const MHDFloat lb = -1.0; + const MHDFloat ub = 1.0; + SparseSM::Chebyshev::LinearMap::I4 i4(nN, nN, lb, ub); + qi = i4.mat(); + } + } +} + +void setForcing(Test& test, Matrix& forcing, const MHDFloat t) +{ + auto ref = test.getForcing(); + const auto& a = test.timeCoeffs(0); + const auto& b = test.timeCoeffs(1); + + for (int i = 0; i < ref.rows(); i++) + { + forcing.row(i).array() = ref(i, 0) * std::cos(a * t); + forcing.row(i).array() += ref(i, 1) * std::sin(a * t); + forcing.row(i).array() += ref(i, 2) * std::cos(b * t); + forcing.row(i).array() += ref(i, 3) * std::sin(b * t); + } +} + +void checkSolution(Test& test, const Matrix& data, const int tIdx) +{ + auto ref = test.getReference(); + + auto err = (data - ref.col(tIdx)).array().abs().maxCoeff(); + std::cerr << test.tEnd(tIdx) << " " << err << std::endl; +} + +void createOperators(Test& test, std::map& ops) +{ + const auto& nN = test.nN; + + if (test.basisId == Test::BasisId::WORLAND) + { + if (test.equationId == Test::EquationId::DIFFUSION) + { + const auto& l = test.l; + const auto& w_a = Polynomial::Worland::worland_chebyshev_t::ALPHA; + const auto& w_db = Polynomial::Worland::worland_chebyshev_t::DBETA; + // Compute model's linear operator (without Tau lines) + ops.insert(std::make_pair(ModelOperator::ImplicitLinear::id(), + DecoupledZSparse())); + SparseSM::Worland::I2Lapl opL(nN, nN, w_a, w_db, l); + ops.find(ModelOperator::ImplicitLinear::id())->second.real() = + opL.mat(); + // Compute model's time operator (Mass matrix) + ops.insert( + std::make_pair(ModelOperator::Time::id(), DecoupledZSparse())); + SparseSM::Worland::I2 opT(nN, nN, w_a, w_db, l); + SparseSM::Worland::Id qid(nN, nN, w_a, w_db, l, nN - 1); + ops.find(ModelOperator::Time::id())->second.real() = + opT.mat() - opT.mat() * qid.mat(); + // Compute model's boundary operator + ops.insert( + std::make_pair(ModelOperator::Boundary::id(), DecoupledZSparse())); + SparseMatrix bc(nN, nN); + std::vector> coeffs; + SparseSM::Worland::Boundary::Value bcValue(w_a, w_db, l); + auto bcVal = bcValue.compute(nN - 1).cast(); + for (int i = 0; i < nN; i++) + { + coeffs.push_back(Eigen::Triplet(0, i, bcVal(i))); + } + bc.setFromTriplets(coeffs.begin(), coeffs.end()); + ops.find(ModelOperator::Boundary::id())->second.real() = bc; + } + else if (test.equationId == Test::EquationId::BIDIFFUSION) + { + const auto& l = test.l; + const auto& w_a = Polynomial::Worland::worland_chebyshev_t::ALPHA; + const auto& w_db = Polynomial::Worland::worland_chebyshev_t::DBETA; + // Compute model's linear operator (without Tau lines) + ops.insert(std::make_pair(ModelOperator::ImplicitLinear::id(), + DecoupledZSparse())); + SparseSM::Worland::I4Lapl2 opL(nN, nN, w_a, w_db, l); + ops.find(ModelOperator::ImplicitLinear::id())->second.real() = + opL.mat(); + // Compute model's time operator (Mass matrix) + ops.insert( + std::make_pair(ModelOperator::Time::id(), DecoupledZSparse())); + SparseSM::Worland::I4Lapl opT(nN, nN, w_a, w_db, l); + ops.find(ModelOperator::Time::id())->second.real() = opT.mat(); + // Compute model's boundary operator + ops.insert( + std::make_pair(ModelOperator::Boundary::id(), DecoupledZSparse())); + SparseMatrix bc(nN, nN); + std::vector> coeffs; + SparseSM::Worland::Boundary::Value bcValue(w_a, w_db, l); + auto bcVal = bcValue.compute(nN - 1).cast(); + SparseSM::Worland::Boundary::D2 bcDiff2(w_a, w_db, l); + auto bcD2 = bcDiff2.compute(nN - 1).cast(); + for (int i = 0; i < nN; i++) + { + coeffs.push_back(Eigen::Triplet(0, i, bcVal(i))); + coeffs.push_back(Eigen::Triplet(1, i, bcD2(i))); + } + bc.setFromTriplets(coeffs.begin(), coeffs.end()); + ops.find(ModelOperator::Boundary::id())->second.real() = bc; + } + } + else if (test.basisId == Test::BasisId::CHEBYSHEV) + { + const MHDFloat lb = -1.0; + const MHDFloat ub = 1.0; + + // Compute model's linear operator (without Tau lines) + ops.insert(std::make_pair(ModelOperator::ImplicitLinear::id(), + DecoupledZSparse())); + SparseMatrix opL(nN, nN); + std::vector> coeffs; + for (int i = 2; i < nN; i++) + { + coeffs.push_back(Eigen::Triplet(i, i, 1.0)); + } + opL.setFromTriplets(coeffs.begin(), coeffs.end()); + ops.find(ModelOperator::ImplicitLinear::id())->second.real() = opL; + // Compute model's time operator (Mass matrix) + ops.insert(std::make_pair(ModelOperator::Time::id(), DecoupledZSparse())); + SparseSM::Chebyshev::LinearMap::I2 opT(nN, nN, lb, ub); + ops.find(ModelOperator::Time::id())->second.real() = opT.mat(); + // Compute model's boundary operator + ops.insert( + std::make_pair(ModelOperator::Boundary::id(), DecoupledZSparse())); + SparseMatrix bc(nN, nN); + coeffs.clear(); + coeffs.push_back(Eigen::Triplet(0, 0, 1.0)); + coeffs.push_back(Eigen::Triplet(1, 0, 1.0)); + for (int i = 1; i < nN; i++) + { + coeffs.push_back(Eigen::Triplet(0, i, 2.0)); + coeffs.push_back( + Eigen::Triplet(1, i, std::pow(-1.0, i) * 2.0)); + } + bc.setFromTriplets(coeffs.begin(), coeffs.end()); + ops.find(ModelOperator::Boundary::id())->second.real() = bc; + } +} + +void initSolution(Test& test, DecoupledZMatrix& sol) +{ + auto ref = test.getInitial(); + sol.real() = ref; + sol.imag() = ref; +} + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace TestSuite +} // namespace QuICC diff --git a/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestHelper.hpp b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestHelper.hpp new file mode 100644 index 00000000..6edbc6e4 --- /dev/null +++ b/Components/TimestepPcViews/TestSuite/Timestep/PredictorCorrector/TestHelper.hpp @@ -0,0 +1,257 @@ +/** + * @file TestHelper.hpp + * @brief Helper functions to setup Timestep tests + */ + +#ifndef QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TESTHELPER_HPP +#define QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TESTHELPER_HPP + +// System includes +// +#include + +// Project includes +// +#include "QuICC/Enums/FieldIds.hpp" +#include "QuICC/PhysicalNames/Temperature.hpp" +#include "QuICC/SolveTiming/Prognostic.hpp" +#include "QuICC/Solver/SparseSolver.hpp" +#include "TestSuite/Timestep/PredictorCorrector/Test.hpp" +#include "Types/Typedefs.hpp" + +namespace QuICC { + +namespace TestSuite { + +namespace Timestep { + +namespace PredictorCorrector { + +void createOperators(Test& test, std::map& ops); +void initQuasiInverse(Test& test, SparseMatrix& qi); +void setForcing(Test& test, Matrix& forcing, const MHDFloat t); +void initSolution(Test& test, DecoupledZMatrix& sol); +void checkSolution(Test& test, const Matrix& data, const int tIdx); + +template +std::shared_ptr> +createTimestepper(Test& test) +{ + typedef typename TScheme::template ImplementationType + Timestepper; + + const auto& start = test.start; + auto spStepper = + std::make_shared(start, SolveTiming::Prognostic::id()); + + auto spScheme = std::make_shared(); + spStepper->setScheme(spScheme); + + return spStepper; +} + +template +void setupTimestepper(Test& test, + std::shared_ptr> + spStepper) +{ + const int nSystems = 1; + const auto& nN = test.nN; + const auto& cols = test.cols; + const auto& startRow = test.startRow; + const auto& dt = test.dt; + const int idx = 0; + + spStepper->addStorage(nN, cols); + auto fieldId = std::make_pair(PhysicalNames::Temperature::id(), + FieldComponents::Spectral::SCALAR); + spStepper->addInformation(fieldId, 0, startRow); + spStepper->initStartRow(); + + spStepper->initMatrices(nSystems); + + // Build operators + std::map ops; + createOperators(test, ops); + spStepper->buildOperators(idx, ops, dt, nN); + + // Solver is initialized + spStepper->setInitialized(); + spStepper->initSolver(); + + // Create initial state + initSolution(test, spStepper->rSolution(idx)); + spStepper->initSolutions(); +} + +template +void singleStepForward(Test& test, + std::shared_ptr> + spStepper) +{ + const auto& nN = test.nN; + const auto& cols = test.cols; + const auto& tEnd = test.tEnd; + const int idx = 0; + + // Initialize quasi-inverse operators + SparseMatrix qi; + initQuasiInverse(test, qi); + + std::pair status = std::make_pair(false, -1.0); + + // Test single steps + for (int tIdx = 0; tIdx < tEnd.size(); tIdx++) + { + MHDFloat t = test.t0; + auto t_ = t; + // Update timestep + const auto& dt = tEnd(tIdx); + spStepper->updateTimeMatrix(dt); + spStepper->updateSolver(); + + // Reset initial state + initSolution(test, spStepper->rSolution(idx)); + spStepper->initSolutions(); + + do + { + // Set RHS + Matrix forcing(nN, cols); + setForcing(test, forcing, t_); + spStepper->rRHSData(idx).real() = qi * forcing; + spStepper->rRHSData(idx).imag() = qi * forcing; + + bool solving = false; + do + { + // Prepare solve of linear system + bool needSolve = spStepper->preSolve(); + + if (needSolve) + { + // Solve linear system + spStepper->solve(); + + // Work on fields after solve + solving = spStepper->postSolve(); + } + else + { + solving = false; + } + + } while (solving); + + status.first = spStepper->finished(); + + if (status.first) + { + status.second = std::max(status.second, spStepper->error()); + t += dt; + t_ = t; + } + else + { + t_ = t + spStepper->stepFraction() * dt; + } + } while (!status.first); + + // Check solution + checkSolution(test, spStepper->solution(idx).real(), tIdx); + } +} + +template +void multiStepForward(Test& test, + std::shared_ptr> + spStepper, + const int steps) +{ + const auto& nN = test.nN; + const auto& cols = test.cols; + const auto& tEnd = test.tEnd; + const int idx = 0; + + // Initialize quasi-inverse operators + SparseMatrix qi; + initQuasiInverse(test, qi); + + std::pair status = std::make_pair(false, -1.0); + + // Test single steps + for (int tIdx = 0; tIdx < tEnd.size(); tIdx++) + { + MHDFloat t = test.t0; + auto t_ = t; + // Update timestep + const auto& dt = tEnd(tIdx) / static_cast(steps); + spStepper->updateTimeMatrix(dt); + spStepper->updateSolver(); + + // Reset initial state + initSolution(test, spStepper->rSolution(idx)); + spStepper->initSolutions(); + + do + { + do + { + // Set RHS + Matrix forcing(nN, cols); + setForcing(test, forcing, t_); + spStepper->rRHSData(idx).real() = qi * forcing; + spStepper->rRHSData(idx).imag() = qi * forcing; + + bool solving = false; + do + { + // Prepare solve of linear system + bool needSolve = spStepper->preSolve(); + + if (needSolve) + { + // Solve linear system + spStepper->solve(); + + // Work on fields after solve + solving = spStepper->postSolve(); + } + else + { + solving = false; + } + + } while (solving); + + status.first = spStepper->finished(); + + if (status.first) + { + status.second = std::max(status.second, spStepper->error()); + t += dt; + t_ = t; + } + else + { + t_ = t + spStepper->stepFraction() * dt; + } + } while (!status.first); + } while (t < tEnd(tIdx)); + + // Check solution + checkSolution(test, spStepper->solution(idx).real(), tIdx); + } +} + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace TestSuite +} // namespace QuICC + +#endif // QUICC_TESTSUITE_TIMESTEP_PREDICTORCORRECTOR_TESTHELPER_HPP diff --git a/Components/TimestepPcViews/Timestep/CMakeLists.txt b/Components/TimestepPcViews/Timestep/CMakeLists.txt new file mode 100644 index 00000000..4a89b581 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(PredictorCorrector) diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/CMakeLists.txt b/Components/TimestepPcViews/Timestep/PredictorCorrector/CMakeLists.txt new file mode 100644 index 00000000..32cf1549 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/CMakeLists.txt @@ -0,0 +1,5 @@ +quicc_target_sources(${QUICC_CURRENT_COMPONENT_LIB} ${QUICC_CMAKE_SRC_VISIBILITY} + FactoryViews.cpp + ) + +add_subdirectory(Views) diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/FactoryViews.cpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/FactoryViews.cpp new file mode 100644 index 00000000..c24c34b3 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/FactoryViews.cpp @@ -0,0 +1,41 @@ +/** + * @file Factory.cpp + * @brief Implementation of factory for Predictor-Corrector schemes + */ + +// System includes +// + +// Project includes +// +#include "QuICC/Timestep/Id/ImexPc2b.hpp" +#include "Timestep/PredictorCorrector/FactoryViews.hpp" +#include "Timestep/PredictorCorrector/ImExPC2.hpp" +#include "Timestep/PredictorCorrector/InterfaceViews.hpp" + +namespace QuICC { + +namespace Timestep { + +namespace PredictorCorrector { + +std::shared_ptr makeInterfaceViews(const std::size_t schemeId, + const MHDFloat time, const Matrix& cfl, const MHDFloat maxError, + const Timestep::Interface::ScalarEquation_range& scalEq, + const Timestep::Interface::VectorEquation_range& vectEq, + Pseudospectral::Coordinator& pseudo) +{ + std::shared_ptr iface; + + if (schemeId == Id::ImexPc2b::id()) + { + iface = std::make_shared>(time, cfl, maxError, scalEq, + vectEq, pseudo); + } + + return iface; +} + +} // namespace PredictorCorrector +} // namespace TimestepViews +} // namespace QuICC diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/FactoryViews.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/FactoryViews.hpp new file mode 100644 index 00000000..024edb93 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/FactoryViews.hpp @@ -0,0 +1,34 @@ +/** + * @file Factory.hpp + * @brief Predictor-Corrector timestepper factory + */ + +#ifndef QUICC_TIMESTEP_PREDICTORCORRECTOR_FACTORYVIEWS_HPP +#define QUICC_TIMESTEP_PREDICTORCORRECTOR_FACTORYVIEWS_HPP + +// System includes +// +#include + +// Project includes +// +#include "QuICC/Timestep/Interface.hpp" + +namespace QuICC { + +namespace Timestep { + +/// @brief This namespace contains predictor-corrector timestepping schemes +namespace PredictorCorrector { + +std::shared_ptr makeInterfaceViews(const std::size_t schemeId, + const MHDFloat time, const Matrix& cfl, const MHDFloat maxError, + const Timestep::Interface::ScalarEquation_range& scalEq, + const Timestep::Interface::VectorEquation_range& vectEq, + Pseudospectral::Coordinator& pseudo); + +} +} // namespace Timestep +} // namespace QuICC + +#endif // QUICC_TIMESTEP_PREDICTORCORRECTOR_FACTORYVIEWS_HPP diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/InterfaceViews.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/InterfaceViews.hpp new file mode 100644 index 00000000..f0a9c6e9 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/InterfaceViews.hpp @@ -0,0 +1,932 @@ +/** + * @file Interface.hpp + * @brief Implementation of an interface to the ImEx predictor-corrector schemes + */ + +#ifndef QUICC_TIMESTEP_PREDICTORCORRECTOR_INTERFACEVIEWS_HPP +#define QUICC_TIMESTEP_PREDICTORCORRECTOR_INTERFACEVIEWS_HPP + +// System includes +// +#include +#include + +// Project includes +// +#include "Memory/MemoryResource.hpp" +#include "Profiler/Interface.hpp" +#include "QuICC/Debug/DebuggerMacro.h" +#include "QuICC/Equations/CouplingInformation.hpp" +#include "QuICC/Pseudospectral/Coordinator.hpp" +#include "QuICC/SolveTiming/Prognostic.hpp" +#include "QuICC/Timestep/Constants.hpp" +#include "QuICC/Timestep/IScheme.hpp" +#include "QuICC/Timestep/Interface.hpp" +#include "Timestep/PredictorCorrector/Views/TimestepperCoordinator.hpp" +#include "Timestep/PredictorCorrector/Views/ImExPCTimestepper.hpp" +#include "Timestep/PredictorCorrector/Views/details/TimesteppperTools.hpp" +#include "QuICC/Tools/Formatter.hpp" +#include "View/ViewDense.hpp" +#include "Memory/Memory.hpp" +#include "Memory/Cpu/NewDelete.hpp" + +namespace QuICC { + +namespace Timestep { + +namespace PredictorCorrector { + +/** + * @brief Implementation of general timestepper structure + */ +template class InterfaceViews : public Timestep::Interface +{ +public: + /// Typedef for solver implementation + template + using SolverImplementationType = + Views::ImExPCTimestepper; + + /// Typedef for parent coordinator + typedef Views::TimestepperCoordinator + SolverCoordinator; + + /** + * @brief Constructor + * + * @param time Initial time value + * @param cfl Initial CFL timestep + * @param error Max error allowed during timestep + * @param scalEq Shared scalar equations + * @param vectEq Shared vector equations + */ + InterfaceViews(const MHDFloat time, const Matrix& cfl, const MHDFloat maxError, + const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq, + Pseudospectral::Coordinator& pseudo); + + /** + * @brief Destructor + */ + virtual ~InterfaceViews() = default; + + /** + * @brief Timestep is finished? + */ + bool finishedStep() const final; + + /** + * @brief Update equation explicit linear input to solver + * + * @param scalEq Scalar equations + * @param vectEq Vector equations + */ + void getExplicitInput(const std::size_t opId, + const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq, + const ScalarVariable_map& scalVar, + const VectorVariable_map& vectVar) final; + + /** + * @brief Tune adaptive timestepper + * + * \mhdBug Not fully implemented + */ + void tuneAdaptive(const MHDFloat time) final; + + /** + * @brief Adapt the timestep used + * + * @param cfl CFL conditions + * @param scalEq Shared scalar equations + * @param vectEq Shared vector equations + */ + void adaptTimestep(const Matrix& cfl, const ScalarEquation_range& scalEq, + const VectorEquation_range& vectEq) final; + + /** + * @brief Compute (partial) forward step + * + * @param scalEq Shared scalar equations + * @param vectEq Shared vector equations + * @param scalVar Shared scalar variables + * @param vectVar Shared vector variables + */ + void stepForward(const ScalarEquation_range& scalEq, + const VectorEquation_range& vectEq, const ScalarVariable_map& scalVar, + const VectorVariable_map& vectVar) final; + + /** + * @brief Print timestepper information to stream + * + * @param stream Output stream + */ + void printInfo(std::ostream& stream) final; + +protected: + /** + * @brief Initialize solution + * + * @param scalEq Scalar equations + * @param vectEq Vector equations + */ + void initSolution(const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq); + + /** + * @brief Update equation input to solver + * + * @param scalEq Scalar equations + * @param vectEq Vector equations + */ + void getInput(const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq); + + /** + * @brief Transfer solution from solver + * + * @param scalEq Scalar equations + * @param vectEq Vector equations + */ + void transferOutput(const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq); + + /** + * @brief Translate equations to timestepper info + * + * @param infos Vector of information + * @param scalEq Scalar equations + * @param vectEq Vector equations + */ + void translate(std::vector& infos, + const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq); + +private: + /** + * @brief Interface to timestepping scheme + */ + SharedIScheme mspScheme; + + /** + * @brief Interface to timestepping scheme + */ + SolverCoordinator mSolverCoord; + + /** + * @brief + */ + std::shared_ptr _mem; +}; + +/** + * @brief Compute flat timestepper index + */ +std::size_t stepperIndex(const std::size_t solverIndex, const std::size_t sysIndex); + +Views::TimestepperInfo createInfo(const Equations::CouplingInformation& cinfo); + +inline Views::TimestepperInfo createInfo(const Equations::CouplingInformation& cinfo, const std::size_t idx) +{ + Views::TimestepperInfo info; + info.isComplex = cinfo.isComplex(); + info.fieldIndex = cinfo.fieldIndex(); + info.solverIndex = stepperIndex(cinfo.solverIndex(), idx); + info.rows = cinfo.systemN(idx); + info.cols = cinfo.rhsCols(idx); + info.blockN = cinfo.galerkinN(idx); + + return info; +} + +/** + * @brief Wrapper to build timestepping matrices + */ +void buildTimestepMatrixWrapper(std::map& ops, Equations::SharedIEquation spEq, FieldComponents::Spectral::Id comp, const int idx); + +template +InterfaceViews::InterfaceViews(const MHDFloat time, const Matrix& cfl, + const MHDFloat maxError, const ScalarEquation_range& scalEq, + const VectorEquation_range& vectEq, Pseudospectral::Coordinator& pseudo) : + Timestep::Interface(time, cfl, maxError, scalEq, vectEq, pseudo) +{ + _mem = std::make_shared(); + + std::cerr << "##############################################" << std::endl; + std::cerr << "##############################################" << std::endl; + std::cerr << "##############################################" << std::endl; + std::cerr << "INITIALIZING NEW TIMESTEPPER INFRASTRUCTURE" << std::endl; + std::cerr << "##############################################" << std::endl; + std::cerr << "##############################################" << std::endl; + std::cerr << "##############################################" << std::endl; + + // Create Timestepper scheme + std::shared_ptr spScheme = std::make_shared(); + + // Use embedded scheme to compute error + if (maxError > 0.0) + { + spScheme->enableEmbedded(); + this->mMaxError = maxError; + } + this->mspScheme = spScheme; + + std::vector infos; + this->translate(infos, scalEq, vectEq); + + this->mSolverCoord.init(this->timestep(), infos, spScheme); + + this->initSolution(scalEq, vectEq); +} + +// \todo convert to free function +template +void InterfaceViews::translate(std::vector& infos, const ScalarEquation_range& scalEq, + const VectorEquation_range& vectEq) +{ + auto addInfo = [](auto& infos, auto&& eq_range) + { + SpectralFieldId myId; + + // Loop over range + for(auto& eqIt: make_range(eq_range)) + { + for(auto& compIt: make_range(eqIt->spectralRange())) + { + // Get field identity + myId = std::make_pair(eqIt->name(), compIt); + + const auto& cinfo = eqIt->couplingInfo(myId.second); + if(eqIt->solveTiming() == SolveTiming::Prognostic::id()) + { + DebuggerMacro_msg("Creating timesteppers for " + PhysicalNames::Coordinator::tag(eqIt->name()) + "(" + Tools::IdToHuman::toString(static_cast(compIt)) + ")", 2); + + for(std::size_t i = cinfo.fieldStart(); i < static_cast(cinfo.nSystems()); i++) + { + auto info = createInfo(cinfo, i); + + // Set operators + buildTimestepMatrixWrapper(info.ops, eqIt, myId.second, i); + + infos.push_back(info); + } + } + } + } + }; + + addInfo(infos, scalEq); + addInfo(infos, vectEq); +} + +inline std::size_t stepperIndex(const std::size_t solverIndex, const std::size_t sysIndex) +{ + assert(sysIndex < 10000); + + return 10000*solverIndex + sysIndex; +} + +template +void InterfaceViews::tuneAdaptive(const MHDFloat time) +{ + this->mStepTime = time; +} + +template bool InterfaceViews::finishedStep() const +{ + return this->mSolverCoord.finishedStep(); +} + +template +void InterfaceViews::initSolution(const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq) +{ + auto processSolution = [this](auto&& eq_range) + { + // Storage for information and identity + SpectralFieldId myId; + + // Loop over all scalar equations + for(auto& eqIt: make_range(eq_range)) + { + for(auto& compIt: make_range(eqIt->spectralRange())) + { + // Get field identity + myId = std::make_pair(eqIt->name(), compIt); + + const auto& cinfo = eqIt->couplingInfo(myId.second); + + if(eqIt->solveTiming() == SolveTiming::Prognostic::id()) + { + for(std::size_t i = cinfo.fieldStart(); i < static_cast(cinfo.nSystems()); i++) + { + auto info = createInfo(cinfo, i); + + DecoupledZMatrix tmp(cinfo.galerkinN(i), cinfo.rhsCols(i)); + tmp.setZero(); + + std::visit( + [&](auto&& p) + { + Equations::copyUnknown(*eqIt, p->dom(0).perturbation(), myId.second, tmp, i, 0, true, true); + }, eqIt->spUnknown()); + + std::uint32_t mem_rows = static_cast(cinfo.galerkinN(i)); + std::uint32_t mem_cols = static_cast(cinfo.rhsCols(i)); + Memory::MemBlock data(mem_rows*mem_cols, this->_mem.get()); + using dense2D = View::DimLevelType; + std::array dimensions {mem_rows, mem_cols}; + View::View> tmpView(data, dimensions); + Views::details::computeSet(tmpView, tmp, 0); + + this->mSolverCoord.updateSolution(info, tmpView); + } + } + } + } + }; + + processSolution(scalEq); + processSolution(vectEq); +} + +template +void InterfaceViews::getExplicitInput(const std::size_t opId, + const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq, + const ScalarVariable_map& scalVar, const VectorVariable_map& vectVar) +{ + Profiler::RegionFixture<2> fix("Timestep-explicitInput"); + + auto processInput = [this](auto&& eq_range, const std::size_t opId, + const ScalarVariable_map& scalVar, const VectorVariable_map& vectVar) + { + // Storage for information and identity + SpectralFieldId myId; + + // Loop over all scalar equations + for(auto& eqIt: make_range(eq_range)) + { + for(auto& compIt: make_range(eqIt->spectralRange())) + { + // Get field identity + myId = std::make_pair(eqIt->name(), compIt); + + const auto& cinfo = eqIt->couplingInfo(myId.second); + + if(eqIt->solveTiming() == SolveTiming::Prognostic::id()) + { + DebuggerMacro_msg("Get explicit timestepper input for " + PhysicalNames::Coordinator::tag(myId.first) + "(" + Tools::IdToHuman::toString(static_cast(myId.second)) + ")", 6); + + // Build range of operator + auto r = make_range(cinfo.explicitRange(opId)); + +#ifdef QUICC_DEBUG + if(r.size() == 0) + { + DebuggerMacro_msg("(Nothing)", 7); + } +#endif // QUICC_DEBUG + + if(r.size() > 0) + { + // Get timestep input + for(std::size_t i = cinfo.fieldStart(); i < static_cast(cinfo.nSystems()); i++) + { + auto info = createInfo(cinfo, i); + + // Copy field values into timestepper input + DecoupledZMatrix tmp(cinfo.galerkinN(i), cinfo.rhsCols(i)); + tmp.setZero(); + + // Loop over explicit fields + for(auto& fIt: r) + { + DebuggerMacro_msg("Add " + ModelOperator::Coordinator::tag(opId) + " term from " + PhysicalNames::Coordinator::tag(fIt.first) + "(" + Tools::IdToHuman::toString(static_cast(fIt.second)) + ")", 7); + + // Get explicit input + if(fIt.second == FieldComponents::Spectral::SCALAR) + { + std::visit( + [&](auto&& p) + { + Equations::addExplicitTerm(*eqIt, opId, myId.second, tmp, 0, fIt, p->dom(0).perturbation(), i); + }, scalVar.find(fIt.first)->second); + } else + { + std::visit( + [&](auto&& p) + { + Equations::addExplicitTerm(*eqIt, opId, myId.second, tmp, 0, fIt, p->dom(0).perturbation().comp(fIt.second), i); + }, vectVar.find(fIt.first)->second); + } + } + + std::uint32_t mem_rows = static_cast(cinfo.galerkinN(i)); + std::uint32_t mem_cols = static_cast(cinfo.rhsCols(i)); + Memory::MemBlock data(mem_rows*mem_cols, this->_mem.get()); + using dense2D = View::DimLevelType; + std::array dimensions {mem_rows, mem_cols}; + View::View> tmpView(data, dimensions); + Views::details::computeSet(tmpView, tmp, 0); + + this->mSolverCoord.updateRhs(info, tmpView); + } + } + } + } + } + }; + + processInput(scalEq, opId, scalVar, vectVar); + processInput(vectEq, opId, scalVar, vectVar); +} + +template +void InterfaceViews::getInput(const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq) +{ + Profiler::RegionFixture<2> fix("Timestep-input"); + + auto processInput = [this](auto&& eq_range) + { + // Storage for information and identity + SpectralFieldId myId; + + // Loop over all scalar equations + for(auto& eqIt: make_range(eq_range)) + { + for(auto& compIt: make_range(eqIt->spectralRange())) + { + // Get field identity + myId = std::make_pair(eqIt->name(), compIt); + + // Apply constraint on solution + eqIt->applyConstraint(myId.second, SolveTiming::Before::id()); + + const auto& cinfo = eqIt->couplingInfo(myId.second); + + if(eqIt->solveTiming() == SolveTiming::Prognostic::id()) + { + DebuggerMacro_msg("Get timestepper input for " + PhysicalNames::Coordinator::tag(myId.first) + "(" + Tools::IdToHuman::toString(static_cast(myId.second)) + ")", 6); + + // Get timestep input + for(std::size_t i = cinfo.fieldStart(); i < static_cast(cinfo.nSystems()); i++) + { + auto info = createInfo(cinfo, i); + + // Copy field values into timestepper input + DecoupledZMatrix tmp(cinfo.galerkinN(i), cinfo.rhsCols(i)); + tmp.setZero(); + Equations::copyNonlinear(*eqIt, myId.second, tmp, i, 0); + + // Add source term + std::visit( + [&](auto&& p) + { + Equations::addSource(*eqIt, p->dom(0).perturbation(), myId.second, tmp, i, 0); + }, eqIt->spUnknown()); + + std::uint32_t mem_rows = static_cast(cinfo.galerkinN(i)); + std::uint32_t mem_cols = static_cast(cinfo.rhsCols(i)); + Memory::MemBlock data(mem_rows*mem_cols, this->_mem.get()); + using dense2D = View::DimLevelType; + std::array dimensions {mem_rows, mem_cols}; + View::View> tmpView(data, dimensions); + Views::details::computeSet(tmpView, tmp, 0); + + this->mSolverCoord.updateRhs(info, tmpView); + + // If required set inhomogenous boundary condition value + if(cinfo.hasBoundaryValue()) + { + // Set boundary value + std::visit( + [&](auto&& p) + { + throw std::logic_error("NOT YET IMPLEMENTED"); + // Equations::setBoundaryValue(*spEq, p->dom(0).perturbation(), id.second, (*solveIt)->rInhomogeneous(i), i, (*solveIt)->startRow(id,i)); + }, eqIt->spUnknown()); + // this->mSolverCoord.updateInhomogeneous(info); + } + } + } + } + } + }; + + processInput(scalEq); + processInput(vectEq); +} + +template +void InterfaceViews::transferOutput(const ScalarEquation_range& scalEq, const VectorEquation_range& vectEq) +{ + Profiler::RegionFixture<2> fix("Timestep-output"); + + auto processOutput = [this](auto&& eq_range) + { + // Storage for information and identity + SpectralFieldId myId; + + // Loop over all scalar equations + for(auto& eqIt: make_range(eq_range)) + { + for(auto& compIt: make_range(eqIt->spectralRange())) + { + // Get field identity + myId = std::make_pair(eqIt->name(), compIt); + + const auto& cinfo = eqIt->couplingInfo(myId.second); + + if(eqIt->solveTiming() == SolveTiming::Prognostic::id()) + { + DebuggerMacro_msg("Get timestepper solution for " + PhysicalNames::Coordinator::tag(myId.first) + "(" + Tools::IdToHuman::toString(static_cast(myId.second)) + ")", 6); + + // return zero + for(std::size_t i = 0; i < static_cast(cinfo.fieldStart()); i++) + { + DecoupledZMatrix tmp(cinfo.galerkinN(i), cinfo.rhsCols(i)); + tmp.setZero(); + + eqIt->storeSolution(myId.second, tmp, i, 0); + } + + for(std::size_t i = cinfo.fieldStart(); i < static_cast(cinfo.nSystems()); i++) + { + auto info = createInfo(cinfo, i); + + std::uint32_t mem_rows = static_cast(cinfo.galerkinN(i)); + std::uint32_t mem_cols = static_cast(cinfo.rhsCols(i)); + Memory::MemBlock data(mem_rows*mem_cols, this->_mem.get()); + using dense2D = View::DimLevelType; + std::array dimensions {mem_rows, mem_cols}; + View::View> tmpView(data, dimensions); + this->mSolverCoord.getSolution(tmpView, info); + + DecoupledZMatrix tmp(cinfo.galerkinN(i), cinfo.rhsCols(i)); + Views::details::computeSet(tmp, tmpView, 0); + + eqIt->storeSolution(myId.second, tmp, i, 0); + } + + // Apply constraint on solution + auto changedSolution = eqIt->applyConstraint(myId.second, SolveTiming::After::id()); + + // Update timestepper solver solution if constraint modified it + if(changedSolution) + { + for(std::size_t i = cinfo.fieldStart(); i < static_cast(cinfo.nSystems()); i++) + { + auto info = createInfo(cinfo, i); + + DecoupledZMatrix tmp(cinfo.galerkinN(i), cinfo.rhsCols(i)); + tmp.setZero(); + + std::visit( + [&](auto&& p) + { + Equations::copyUnknown(*eqIt, p->dom(0).perturbation(), myId.second, tmp, i, 0, true, true); + }, eqIt->spUnknown()); + + std::uint32_t mem_rows = static_cast(cinfo.galerkinN(i)); + std::uint32_t mem_cols = static_cast(cinfo.rhsCols(i)); + Memory::MemBlock data(mem_rows*mem_cols, this->_mem.get()); + using dense2D = View::DimLevelType; + std::array dimensions {mem_rows, mem_cols}; + View::View> tmpView(data, dimensions); + Views::details::computeSet(tmpView, tmp, 0); + + this->mSolverCoord.updateSolution(info, tmpView); + } + } + } + } + } + }; + + processOutput(scalEq); + processOutput(vectEq); +} + +template +void InterfaceViews::adaptTimestep(const Matrix& cfl, + const ScalarEquation_range&, const VectorEquation_range&) +{ + // Store old timestep + this->mOldDt = this->timestep(); + + // Update CFL information + this->mDt.block(0, 1, this->mDt.rows(), cfl.cols() - 1) = + cfl.rightCols(cfl.cols() - 1); + + // New computed CFL + MHDFloat compCfl = cfl(0, 0); + + // Check if CFL allows for a larger timestep + MHDFloat newCflDt = 0.0; + if (compCfl > this->mcUpWindow * this->timestep()) + { + if (this->mCnstSteps >= this->mcMinCnst) + { + // Set new timestep + newCflDt = std::min(compCfl, this->mcMaxJump * this->timestep()); + } + else + { + // Reuse same timestep + newCflDt = this->timestep(); + } + + // Check if CFL is below minimal timestep or downard jump is large + } + else if (compCfl < this->mcMinDt || + compCfl < this->timestep() / this->mcMaxJump) + { + // Signal simulation abort + newCflDt = -compCfl; + + // Check if CFL requires a lower timestep + } + else if (compCfl < this->timestep() * (2.0 - this->mcUpWindow)) + { + // Set new timestep + newCflDt = compCfl; + } + else + { + newCflDt = this->timestep(); + } + + // Get timestepper error (if applicable) + MHDFloat error = this->mSolverCoord.error(); + +// Gather error across processes +#ifdef QUICC_MPI + if (error > 0.0) + { + MPI_Allreduce(MPI_IN_PLACE, &error, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); + } +#endif // QUICC_MPI + + // No error control and no CFL condition + MHDFloat newErrorDt = 0.0; + + // Use what ever condition is used by CFL + if (error < 0) + { + newErrorDt = -1.0; + + // Error is too large, reduce timestep + } + else if (error > this->mMaxError) + { + newErrorDt = + this->timestep() * + std::pow(this->mMaxError / error, 1. / this->mspScheme->order()) / + this->mcUpWindow; + + // Error is small, increase timestep + } + else if (error < this->mMaxError / (this->mcMaxJump * 0.9) && + this->mCnstSteps >= this->mcMinCnst) + { + newErrorDt = + std::min(this->timestep() * std::pow(this->mMaxError / error, + 1. / this->mspScheme->order()), + this->timestep() * this->mcMaxJump); + + // Timestep should not be increased + } + else + { + newErrorDt = this->timestep(); + } + + // Update error details + if (this->mMaxError > 0.0) + { + this->mDt(0, this->mDt.cols() - 1) = newErrorDt; + this->mDt(1, this->mDt.cols() - 1) = error; + } + + // CFL condition requested abort! + if (newCflDt < 0.0) + { + this->mDt(0, 0) = newCflDt; + this->mDt(1, 0) = cfl(1, 0); + + // Get minimum between both conditions + } + else if (newCflDt > 0.0 && newErrorDt > 0.0) + { + if (newCflDt < newErrorDt) + { + this->mDt(0, 0) = newCflDt; + this->mDt(1, 0) = cfl(1, 0); + } + else + { + this->mDt(0, 0) = newErrorDt; + this->mDt(1, 0) = ERROR_LOCATION; + } + + // Use CFL condition + } + else if (newCflDt > 0.0) + { + if (this->timestep() != newCflDt) + { + this->mDt(0, 0) = newCflDt; + this->mDt(1, 0) = cfl(1, 0); + } + + // Use error condition + } + else if (newErrorDt > 0.0) + { + this->mDt(0, 0) = newErrorDt; + this->mDt(1, 0) = ERROR_LOCATION; + } + + // + // Update the timestep matrices if necessary + // + if (this->timestep() != this->mOldDt && this->timestep() > 0.0) + { + DebuggerMacro_showValue( + "Updating timestep and matrices with new Dt = ", 0, this->timestep()); + + this->mSolverCoord.updateTimestep(this->timestep()); + + // Update the time dependence in matrices + DebuggerMacro_start("Update matrices", 0); + this->mSolverCoord.updateMatrices(); + DebuggerMacro_stop("Update matrices t = ", 0); + } + else + { + this->mCnstSteps += 1.0; + } + + // Update CFL writer + this->mspIo->setSimTime(this->mTime, this->mDt, this->mCnstSteps); + this->mspIo->write(); + + if (this->timestep() != this->mOldDt && this->timestep() > 0.0) + { + this->mCnstSteps = 0.0; + } +} + +template +void InterfaceViews::stepForward(const ScalarEquation_range& scalEq, + const VectorEquation_range& vectEq, const ScalarVariable_map& scalVar, + const VectorVariable_map& vectVar) +{ + bool isIntegrating = true; + while (isIntegrating) + { + DebuggerMacro_msg("Time integration sub-step", 2); + + this->mpPseudo->evolveUntilPrognostic(this->finishedStep()); + + // Update the equation input to the timestepper + this->getInput(scalEq, vectEq); + + Profiler::RegionStart<2>("Timestep-solve"); + // Solve all the linear systems + this->mSolverCoord.solveSystems(); + Profiler::RegionStop<2>("Timestep-solve"); + + // Transfer timestep output back to equations + this->transferOutput(scalEq, vectEq); + + // Clear the solver RHS + this->mSolverCoord.clearSolvers(); + + // Update current time + this->mTime = + this->mRefTime + this->mSolverCoord.stepFraction() * this->timestep(); + + this->mpPseudo->evolveAfterPrognostic(this->finishedStep()); + + isIntegrating = !this->finishedStep(); + } +} + +template +void InterfaceViews::printInfo(std::ostream& stream) +{ + // Create nice looking ouput header + Tools::Formatter::printNewline(stream); + Tools::Formatter::printLine(stream, '-'); + Tools::Formatter::printCentered(stream, "Timestepper information", '*'); + Tools::Formatter::printLine(stream, '-'); + + std::stringstream oss; + int base = 20; + + // Timestep scheme + oss << "Timestepper: " << this->mspScheme->name() << " (" + << this->mspScheme->order() << ")"; + Tools::Formatter::printCentered(stream, oss.str(), ' ', base); + oss.str(""); + + // General linear solver + oss << "General solver: "; +#if defined QUICC_SPLINALG_MUMPS + oss << "MUMPS"; +#elif defined QUICC_SPLINALG_UMFPACK + oss << "UmfPack"; +#elif defined QUICC_SPLINALG_SPARSELU + oss << "SparseLU"; +#else + oss << "(unknown)"; +#endif // defined QUICC_SPLINALG_MUMPS + + Tools::Formatter::printCentered(stream, oss.str(), ' ', base); + oss.str(""); + + // Triangular linear solver + oss << "Triangular solver: "; +#if defined QUICC_SPTRILINALG_SPARSELU + oss << "SparseLU"; +#elif defined QUICC_SPTRILINALG_MUMPS + oss << "MUMPS"; +#elif defined QUICC_SPTRILINALG_UMFPACK + oss << "UmfPack"; +#else + oss << "(unknown)"; +#endif // defined QUICC_SPTRILINALG_SPARSELU + + Tools::Formatter::printCentered(stream, oss.str(), ' ', base); + oss.str(""); + + // SPD linear solver + oss << "SPD solver: "; +#if defined QUICC_SPSPDLINALG_SIMPLICIALLDLT + oss << "SimplicialLDLT"; +#elif defined QUICC_SPSPDLINALG_SIMPLICIALLLT + oss << "SimplicialLLT"; +#elif defined QUICC_SPSPDLINALG_MUMPS + oss << "MUMPS"; +#elif defined QUICC_SPSPDLINALG_UMFPACK + oss << "UmfPack"; +#elif defined QUICC_SPSPDLINALG_SPARSELU + oss << "SparseLU"; +#else + oss << "(unknown)"; +#endif // defined QUICC_SPSPDLINALG_SIMPLICIALLDLT + + Tools::Formatter::printCentered(stream, oss.str(), ' ', base); + oss.str(""); + + Tools::Formatter::printLine(stream, '*'); + Tools::Formatter::printNewline(stream); +} + +// \todo move to separate file +inline void buildTimestepMatrixWrapper(std::map& ops, Equations::SharedIEquation spEq, FieldComponents::Spectral::Id comp, + const int idx) +{ + bool isSplit = spEq->couplingInfo(comp).isSplitEquation(); + + // Compute model's linear operator (without Tau lines) + ops.insert( + std::make_pair(ModelOperator::ImplicitLinear::id(), DecoupledZSparse())); + spEq->buildModelMatrix(ops.find(ModelOperator::ImplicitLinear::id())->second, + ModelOperator::ImplicitLinear::id(), comp, idx, + ModelOperatorBoundary::SolverNoTau::id()); + // Compute model's time operator (without Tau lines) + ops.insert(std::make_pair(ModelOperator::Time::id(), DecoupledZSparse())); + spEq->buildModelMatrix(ops.find(ModelOperator::Time::id())->second, + ModelOperator::Time::id(), comp, idx, + ModelOperatorBoundary::SolverNoTau::id()); + // Compute model's tau line boundary operator + ops.insert( + std::make_pair(ModelOperator::Boundary::id(), DecoupledZSparse())); + spEq->buildModelMatrix(ops.find(ModelOperator::Boundary::id())->second, + ModelOperator::Boundary::id(), comp, idx, + ModelOperatorBoundary::SolverHasBc::id()); + + // If equation was split into two lower order systems + if (isSplit) + { + // Compute model's split linear operator (without Tau lines) + auto id = ModelOperator::SplitImplicitLinear::id(); + ops.insert(std::make_pair(id, DecoupledZSparse())); + spEq->buildModelMatrix(ops.find(id)->second, id, comp, idx, + ModelOperatorBoundary::SolverNoTau::id()); + + // Compute model's tau line boundary operator for split operator + id = ModelOperator::SplitBoundary::id(); + ops.insert(std::make_pair(id, DecoupledZSparse())); + spEq->buildModelMatrix(ops.find(id)->second, id, comp, idx, + ModelOperatorBoundary::SolverHasBc::id()); + + // Compute model's tau line boundary value for split operator + id = ModelOperator::SplitBoundaryValue::id(); + ops.insert(std::make_pair(id, DecoupledZSparse())); + spEq->buildModelMatrix(ops.find(id)->second, id, comp, idx, + ModelOperatorBoundary::SolverNoTau::id()); + } +} + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace QuICC + +#endif // QUICC_TIMESTEP_PREDICTORCORRECTOR_INTERFACEVIEWS_HPP diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/CMakeLists.txt b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/CMakeLists.txt new file mode 100644 index 00000000..fed7e403 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/CMakeLists.txt @@ -0,0 +1,2 @@ +#quicc_target_sources(${QUICC_CURRENT_COMPONENT_LIB} ${QUICC_CMAKE_SRC_VISIBILITY} +# ) diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ITimestepper.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ITimestepper.hpp new file mode 100644 index 00000000..07829c60 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ITimestepper.hpp @@ -0,0 +1,291 @@ +/** + * @file ITimestepper.hpp + * @brief Implementation of base for the templated (coupled) equation + * timestepper + */ + +#ifndef QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_ITIMESTEPPER_HPP +#define QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_ITIMESTEPPER_HPP + +// System includes +// +#include +#include +#include +#include + +// Project includes +// +#include "QuICC/ModelOperator/Boundary.hpp" +#include "QuICC/ModelOperator/ImplicitLinear.hpp" +#include "QuICC/ModelOperator/SplitBoundary.hpp" +#include "QuICC/ModelOperator/SplitBoundaryValue.hpp" +#include "QuICC/ModelOperator/SplitImplicitLinear.hpp" +#include "QuICC/ModelOperator/Time.hpp" +#include "QuICC/Register/Implicit.hpp" +#include "QuICC/Register/Influence.hpp" +#include "QuICC/Tag/Operator/Rhs.hpp" +#include "Timestep/PredictorCorrector/Views/ITimestepperBase.hpp" +#include "Timestep/PredictorCorrector/Views/details/TimesteppperTools.hpp" +#include "QuICC/Tag/Operator/Influence.hpp" + +namespace QuICC { + +namespace Timestep { + +namespace PredictorCorrector { + +namespace Views { + +/** + * @brief Implementation of a templated (coupled) equation timestepper + */ +template +class ITimestepper + : public ITimestepperBase +{ +public: + /** + * @brief Constructor + */ + ITimestepper(); + + /** + * @brief Destructor + */ + virtual ~ITimestepper() = default; + + /** + * @brief Initialise the solver matrices storage + */ + virtual void initMatrices(); + + /** + * @brief Update the LHS matrix with new timedependence + */ + void updateTimeMatrix(const MHDFloat dt); + + /** + * @brief Finished timestep? + */ + MHDFloat error() const; + + /** + * @brief Finished timestep? + */ + bool finished(); + + /** + * @brief Get current timestep fraction + */ + virtual MHDFloat stepFraction() const = 0; + +protected: + using ITimestepperBase::initMatrices; + + /** + * @brief Number of substeps + */ + virtual int steps() const = 0; + + /** + * @brief + */ + virtual void postSolverUpdate(); + + /** + * @brief Implicit coefficient a for linear operator + * + * A = (T + a L) + */ + virtual MHDFloat aIm(const int step) const = 0; + + /** + * @brief Explicit calculation took place? + */ + bool mHasExplicit; + + /** + * @brief Current substep + */ + int mStep; + + /** + * @brief Current timestep + */ + MHDFloat mDt; + + /** + * @brief Timestep error + */ + MHDFloat mError; + + /** + * @brief ID of the register to use + */ + std::size_t mRegisterId; + + /** + * @brief Mass matrix operator + */ + SparseMatrix mMassMatrix; + +private: +}; + +template +ITimestepper::ITimestepper() : + ITimestepperBase(), + mHasExplicit(true), + mStep(0), + mDt(-1.0), + mError(-1.0), + mRegisterId(Register::Implicit::id()) +{} + +template +MHDFloat ITimestepper::error() const +{ + return this->mError; +} + +template +bool ITimestepper::finished() +{ + return (this->mStep == 0); +} + +template +void ITimestepper::updateTimeMatrix( + const MHDFloat dt) +{ + // Update stored timestep + MHDFloat oldDt = this->mDt; + this->mDt = dt; + + // Get list of different IDs + std::set filter; + for (int step = 0; step < this->steps(); ++step) + { + const auto a = this->aIm(step); + filter.insert(a); + } + + // Loop over all operator IDs + for (auto opIt = this->mOperators.begin(); + opIt != this->mOperators.end(); ++opIt) + { + // Loop of step IDs + for (auto a: filter) + { + // Update is only required if aIm is not zero + if (opIt->second.count(a) > 0 && a != 0.0) + { + // Get the number of nonzero elements in time dependence + size_t nnz = this->linearOperator(Tag::Operator::Rhs::id(), 0).nonZeros(); + + // Update LHS and RHS matrices + for (size_t k = 0; + k < static_cast(this->linearOperator(Tag::Operator::Rhs::id(), 0).outerSize()); + ++k) + { + typename TOperator::InnerIterator lhsIt( + this->linearOperator(opIt->first, a), k); + for (typename TOperator::InnerIterator timeIt( + this->linearOperator(Tag::Operator::Rhs::id(), 0), k); + timeIt; ++timeIt) + { + // Only keep going if nonzero elements are left + if (nnz > 0) + { + assert(lhsIt.col() == timeIt.col()); + assert(lhsIt.row() <= timeIt.row()); + + // LHS matrix might have additional nonzero entries + while (lhsIt.row() < timeIt.row() && lhsIt) + { + ++lhsIt; + } + + // Update LHS matrix + if (timeIt.row() == lhsIt.row()) + { + // Update values + lhsIt.valueRef() += + a * (oldDt - this->mDt) * timeIt.value(); + + // Update nonzero counter + nnz--; + } + + // Update LHS iterators and counters + ++lhsIt; + } + else + { + break; + } + } + } + + // Abort if some nonzero entries where not updated + if (nnz != 0) + { + throw std::logic_error( + "Update of timestepping matrices failed"); + } + } + } + } +} + +template +void ITimestepper::initMatrices() +{ + // Initialise base matrices + for (int i = 0; i < this->steps(); i++) + { + this->initMatrices(Tag::Operator::Lhs::id(), this->aIm(i)); + } + + this->initMatrices(Tag::Operator::Rhs::id(), 0); +} + +template +void ITimestepper::postSolverUpdate() +{ + const auto lhsId = Tag::Operator::Lhs::id(); + const auto opId = Tag::Operator::Influence::id(); + + if (this->mSolver.count(lhsId) > 0 && this->mSolver.count(opId) > 0) + { + // Solver for both stages + auto&& sIt1 = this->mSolver.at(opId).find(0.0)->second; + auto&& sIt2 = this->mSolver.at(lhsId).find(0.0)->second; + + // Compute green's functions + auto&& infKernel = this->reg(Register::Influence::id()); + assert(infKernel.real().rows() == infKernel.imag().rows()); + assert(infKernel.real().cols() == infKernel.imag().cols()); + auto rows = infKernel.real().rows(); + auto cols = infKernel.real().cols() / 3; + TData rhs(rows, cols); + for (int i = 0; i < cols; i++) + { + rhs.real().col(i) = infKernel.real().col(3 * i + 2); + rhs.imag().col(i) = infKernel.imag().col(3 * i + 2); + } + TData sol(rows, cols); + details::solveWrapper(sol, sIt1, rhs); + details::computeMV(rhs, this->mMassMatrix, sol); + details::solveWrapper(sol, sIt2, rhs); + details::computeSetInfluence(infKernel, sol); + } +} + +} // namespace Views +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace QuICC + +#endif // QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_ITIMESTEPPER_HPP diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ITimestepperBase.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ITimestepperBase.hpp new file mode 100644 index 00000000..c1c68185 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ITimestepperBase.hpp @@ -0,0 +1,486 @@ +/** + * @file ITimestepperBase.hpp + * @brief Implementation of a templated (coupled) linear solver structure + */ + +#ifndef QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_ITIMESTEPPERBASE_HPP +#define QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_ITIMESTEPPERBASE_HPP + +// System includes +// +#include +#include + +// Project includes +// +#include "Types/Math.hpp" +#include "QuICC/ModelOperator/ImplicitLinear.hpp" +#include "QuICC/ModelOperator/Boundary.hpp" +#include "QuICC/Framework/MpiFramework.hpp" +#include "QuICC/Solver/SparseSolver.hpp" +#include "Timestep/PredictorCorrector/Views/details/SparseLinearSolverTools.hpp" +#include "Timestep/PredictorCorrector/Views/details/TimesteppperTools.hpp" +#include "QuICC/Register/Solution.hpp" +#include "QuICC/Register/Rhs.hpp" +#include "QuICC/Tag/Operator/Lhs.hpp" + +namespace QuICC { + +namespace Timestep { + +namespace PredictorCorrector { + +namespace Views { + + template class ITimestepper; + + /** + * @brief Implementation of a generic timestepper + */ + template class ITimestepperBase + { + public: + /** + * @brief Constructor + */ + ITimestepperBase(); + + /** + * @brief Destructor + */ + virtual ~ITimestepperBase() = default; + + /** + * @brief Is operator initialized? + */ + bool isInitialized() const; + + /** + * @brief Set operator to initialized + */ + void setInitialized(); + + /** + * @brief Initialise solver + */ + void initSolver(); + + /** + * @brief Update solver + */ + void updateSolver(); + + /** + * @brief Set solver RHS data to zero + */ + void zeroSolver(); + + /** + * @brief Solve linear systems + */ + void solve(); + + /** + * @brief Add RHS and solution data storage + * + * @param rows Number of rows of matrix + * @param cols Number of columns required + */ + virtual void addStorage(const int rows, const int cols) = 0; + + /** + * @brief Init storage for inhomogeneous BC + * + * @param rows Number of rows of matrix + * @param cols Number of columns required + */ + virtual void initInhomogeneous(const int rows, const int cols); + + /** + * @brief Get inhomogeneous boundary condition + */ + const Eigen::SparseMatrix& inhomogeneous() const; + + /** + * @brief Set inhomogeneous boundary condition + */ + Eigen::SparseMatrix& rInhomogeneous(); + + /** + * @brief Add inhomogeneous boundary condition to RHS data + */ + void addInhomogeneous(); + + /** + * @brief Computation error + */ + virtual MHDFloat error() const; + + /** + * @brief Finished computation? + */ + virtual bool finished(); + + protected: + /// Typedef for shared solver + typedef Framework::Selector::SparseSolver SolverType; + + /// Typedef for shared solver + typedef std::shared_ptr SharedSolverType; + + /** + * @brief Has solver matrix? + * + * @param opId Operator ID + */ + bool hasLinearOperator(const std::size_t opId) const; + + /** + * @brief Has solver matrix? + * + * @param opId Operator ID + * @param id ID of matrix + */ + bool hasLinearOperator(const std::size_t opId, const MHDFloat id) const; + + /** + * @brief Set linear operator + * + * @param opId Operator ID + * @param id ID of matrix + */ + TOperator& linearOperator(const std::size_t opId, const MHDFloat id); + + /** + * @brief Get linear operator + * + * @param opId Operator ID + * @param id ID of matrix + */ + const TOperator& linearOperator(const std::size_t opId, const MHDFloat id) const; + + /** + * @brief Get storage register + * + * @param id Register ID + */ + TData& reg(const std::size_t regId); + + /** + * @brief Get storage register + * + * @param id Register ID + */ + const TData& reg(const std::size_t regId) const; + + /** + * @brief Add RHS and solution data storage + * + * @param rows Number of rows of matrix + * @param cols Number of columns required + * @param ids Register IDs + */ + void addRegister(const int rows, const int cols, const std::vector& ids); + + /** + * @brief Initialise the solver matrices storage + * + * @param opId Operator ID + * @param id Id of matrix + */ + void initMatrices(const std::size_t opId, const MHDFloat id); + + /** + * @brief Initialise solver for given ID + * + * @param opId Operator ID + */ + void initSolver(const std::size_t opId); + + /** + * @brief Update solver for given ID + * + * @param opId Operator ID + */ + void updateSolver(const std::size_t opId); + + /** + * @brief Additional updates after solver update + */ + virtual void postSolverUpdate(); + + /** + * @brief Correct solution obtained from linear solver + */ + virtual int correctSolution(const int iteration); + + /** + * @brief Flag for operator initialization + */ + bool mIsInitialized; + + /** + * @brief Current operator ID of the solvers + */ + std::size_t mOpId; + + /** + * @brief Current ID of the solvers + */ + MHDFloat mId; + + private: + /// ITimestepper class is a friend + friend class ITimestepper; + + /** + * @brief Operators used by solver + */ + std::map > mOperators; + + /** + * @brief Create sparse solvers + */ + std::map > mSolver; + + /** + * @brief Storage for field + */ + std::map mStorage; + + /** + * @brief Storage for inhomogeneous boundary conditions + */ + Eigen::SparseMatrix mInhomogeneous; + }; + + template ITimestepperBase::ITimestepperBase() + : mIsInitialized(false), mOpId(Tag::Operator::Lhs::id()), mId(0.0) + { + } + + template MHDFloat ITimestepperBase::error() const + { + return -1.0; + } + + template bool ITimestepperBase::finished() + { + return true; + } + + template TData& ITimestepperBase::reg(const std::size_t id) + { + assert(this->mStorage.count(id) > 0); + + return this->mStorage.at(id); + } + + template const TData& ITimestepperBase::reg(const std::size_t id) const + { + assert(this->mStorage.count(id) > 0); + + return this->mStorage.at(id); + } + + template void ITimestepperBase::addRegister(const int rows, const int cols, const std::vector& ids) + { + // Assert for non zero rows and columns + assert(rows > 0); + assert(cols > 0); + + for(auto id: ids) + { + this->mStorage.emplace(id, TData(rows, cols)); + } + } + + template void ITimestepperBase::addInhomogeneous() + { + if(this->mInhomogeneous.nonZeros() > 0) + { + details::addCorrection(this->reg(Register::Rhs::id()), this->mInhomogeneous); + } + } + + template void ITimestepperBase::solve() + { + int iteration = 0; + + while(iteration >= 0) + { + // Solve other modes + assert(this->mSolver.count(this->mOpId) > 0); + assert(this->mSolver.at(this->mOpId).count(this->mId) > 0); + auto&& spSolver = this->mSolver.at(this->mOpId).find(this->mId)->second; + details::solveWrapper(this->reg(Register::Solution::id()), spSolver, this->reg(Register::Rhs::id())); + + // Stop simulation if solve failed + if(spSolver->info() != Eigen::Success) + { + throw std::logic_error("Sparse direct solve failed!"); + } + + // Callback site for correcting solve solution (for example for influence matrix approach) + iteration = this->correctSolution(iteration); + } + } + + template int ITimestepperBase::correctSolution(const int iteration) + { + return -1; + } + + template void ITimestepperBase::zeroSolver() + { + this->reg(Register::Rhs::id()).setZero(); + } + + template void ITimestepperBase::initSolver() + { + this->initSolver(Tag::Operator::Lhs::id()); + } + + template void ITimestepperBase::initSolver(const std::size_t opId) + { + // Loop over matrices + assert(this->mOperators.count(opId) > 0); + auto&& lhsMatrix = this->mOperators.at(opId); + + // Create operator map + if(this->mSolver.count(opId) == 0) + { + this->mSolver.emplace(opId, std::map()); + } + + // Initialize solvers with matrices + for(auto it = lhsMatrix.begin(); it != lhsMatrix.end(); ++it) + { + auto&& solvers = this->mSolver.at(opId); + + auto spSolver = std::make_shared(); + solvers.emplace(it->first, spSolver); + } + + // Compute pattern and factorisation + this->updateSolver(opId); + } + + template void ITimestepperBase::updateSolver() + { + this->updateSolver(Tag::Operator::Lhs::id()); + } + + template void ITimestepperBase::updateSolver(const std::size_t opId) + { + assert(this->mOperators.count(opId) > 0); + assert(this->mSolver.count(opId) > 0); + + auto&& lhsMatrix = this->mOperators.at(opId); + for(auto it = lhsMatrix.begin(); it != lhsMatrix.end(); ++it) + { + // Compute factorisation + auto&& solver = this->mSolver.at(opId); + typename std::map::iterator sIt = solver.find(it->first); + // Safety assert to make sur matrix is compressed + assert(it->second.isCompressed()); + + sIt->second->compute(it->second); + + // Stop simulation if factorization failed + if(sIt->second->info() != Eigen::Success) + { + throw std::logic_error("Matrix factorization failed!"); + } + } + + // Provide call site for additional updates + this->postSolverUpdate(); + } + + template void ITimestepperBase::postSolverUpdate() + { + // Default implementation does nothing + } + + template void ITimestepperBase::initMatrices(const std::size_t opId, const MHDFloat id) + { + // Create operator ID + if(this->mOperators.count(opId) == 0) + { + this->mOperators.emplace(opId, std::map()); + } + + // Do not reinitialise if work already done by other field + auto&& lhsMatrix = this->mOperators.at(opId); + if(lhsMatrix.count(id) == 0) + { + lhsMatrix.emplace(id, TOperator()); + } + } + + template void ITimestepperBase::initInhomogeneous(const int rows, const int cols) + { + // Assert for non zero rows and columns + assert(rows > 0); + assert(cols > 0); + + // Add storage for inhomogeneous boundary value + this->mInhomogeneous = Eigen::SparseMatrix(rows,cols); + this->mInhomogeneous.setZero(); + } + + template bool ITimestepperBase::hasLinearOperator(const std::size_t opId) const + { + bool res = (this->mOperators.count(opId) > 0); + return res; + } + + template bool ITimestepperBase::hasLinearOperator(const std::size_t opId, const MHDFloat id) const + { + bool res = (this->hasLinearOperator(opId) && this->mOperators.at(opId).count(id) > 0); + return res; + } + + template TOperator& ITimestepperBase::linearOperator(const std::size_t opId, const MHDFloat id) + { + assert(this->mOperators.count(opId) > 0); + assert(this->mOperators.at(opId).count(id) > 0); + + return this->mOperators.at(opId).at(id); + } + + template const TOperator& ITimestepperBase::linearOperator(const std::size_t opId, const MHDFloat id) const + { + assert(this->mOperators.count(opId) > 0); + assert(this->mOperators.at(opId).count(id) > 0); + + return this->mOperators.at(opId).at(id); + } + + template const Eigen::SparseMatrix& ITimestepperBase::inhomogeneous() const + { + return this->mInhomogeneous; + } + + template Eigen::SparseMatrix& ITimestepperBase::rInhomogeneous() + { + return this->mInhomogeneous; + } + + template bool ITimestepperBase::isInitialized() const + { + return this->mIsInitialized; + } + + template void ITimestepperBase::setInitialized() + { + this->mIsInitialized = true; + } + +} // Views +} // PredictorCorrector +} // Timestep +} // QuICC + +#endif // QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_ITIMESTEPPERBASE_HPP diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ImExPCTimestepper.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ImExPCTimestepper.hpp new file mode 100644 index 00000000..d75895e1 --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/ImExPCTimestepper.hpp @@ -0,0 +1,446 @@ +/** + * @file ImExPCTimestepper.hpp + * @brief Implementation of a templated (coupled) equation timestepper for + * Implicit-Explicit Predictor-Corrector schemes. + */ + +#ifndef QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_IMEXPCTIMSTEPPER_HPP +#define QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_IMEXPCTIMSTEPPER_HPP + +// System includes +// + +// Project includes +// +#include "QuICC/Register/Error.hpp" +#include "QuICC/Register/Explicit.hpp" +#include "QuICC/Register/Implicit.hpp" +#include "QuICC/Register/Intermediate.hpp" +#include "QuICC/Register/Solution.hpp" +#include "Timestep//PredictorCorrector/Views/ITimestepper.hpp" +#include "Timestep/PredictorCorrector/IImExPCScheme.hpp" +#include "View/ViewDense.hpp" + +namespace QuICC { + +namespace Timestep { + +namespace PredictorCorrector { + +namespace Views { + +/** + * @brief Implementation of a templated (coupled) equation timestepper for + * Implicit-Explicit Predictor-Corrector schemes + */ +template +class ImExPCTimestepper + : public ITimestepper +{ +public: + /** + * @brief Constructor + */ + ImExPCTimestepper(); + + /** + * @brief Destructor + */ + virtual ~ImExPCTimestepper() = default; + + /** + * @brief Set timestepper scheme + */ + void setScheme(SharedIImExPCScheme spScheme); + + /** + * @brief Number of substeps + */ + int steps() const final; + + /** + * @brief Implicit coefficient a for linear operator + * + * A = (T + a L) + */ + MHDFloat aIm(const int step) const final; + + /** + * + * @brief Get current timestep fraction + */ + virtual MHDFloat stepFraction() const; + + /** + * @brief Prepare fields for implicit solve + */ + bool preSolve(); + + /** + * @brief Work on fields after implicit solve + * + * @param step Current substep + */ + bool postSolve(); + + /** + * @brief Add RHS and solution data storage + * + * @param rows Number of rows of matrix + * @param cols Number of columns required + */ + virtual void addStorage(const int rows, const int cols); + + /** + * @brief Add to RHS + * + * @param rhs RSH values + * @param startRow Start row + */ + void addRhs(const View::View>>& rhs, const std::size_t startRow); + + /** + * @brief Extract solution + * + * @param sol solution values + * @param startRow Start row + */ + void getSolution(View::View>>& sol, const std::size_t startRow); + + /** + * @brief Set initial solution + * + * @param sol initial solution data + * @param startRow Start row + */ + void setSolution(const View::View>>& sol, const std::size_t startRow); + + /** + * @brief Update solver after solution was updated + */ + void updateSolutions(); + + /** + * @brief Build the scheme operators + * + * @param ops Operators for the timestepper + */ + virtual void buildOperators(const std::map& ops, const MHDFloat dt, + const int size); + +protected: + /** + * @brief Timestepping scheme + */ + SharedIImExPCScheme mspScheme; + +private: +}; + +template +ImExPCTimestepper::ImExPCTimestepper() : + ITimestepper() +{} + +template +void ImExPCTimestepper::setScheme( + SharedIImExPCScheme spScheme) +{ + this->mspScheme = spScheme; +} + +template +int ImExPCTimestepper::steps() const +{ + return this->mspScheme->steps(); +} + +template +MHDFloat +ImExPCTimestepper::stepFraction() const +{ + return this->mspScheme->cEx(this->mStep); +} + +template +MHDFloat ImExPCTimestepper::aIm( + const int step) const +{ + return this->mspScheme->aIm(step); +} + +template +void ImExPCTimestepper::updateSolutions() +{ + details::computeSet(this->reg(Register::Intermediate::id()), + this->reg(Register::Solution::id())); + + if (this->mspScheme->useEmbedded()) + { + details::computeSet(this->reg(Register::Error::id()), + this->reg(Register::Solution::id())); + } +} + +template +void ImExPCTimestepper::addStorage( + const int rows, const int cols) +{ + // Assert for non zero rows and columns + assert(rows > 0); + assert(cols > 0); + + // Add additional registers + std::vector ids = { + Register::Solution::id(), + Register::Rhs::id(), + Register::Error::id(), + Register::Explicit::id(), + Register::Intermediate::id() + }; + this->addRegister(rows, cols, ids); + + // Init storage for inhomogeneous boundary value + this->initInhomogeneous(rows, cols); + + // Register for influence kernels + ids = {Register::Influence::id()}; + this->addRegister(1, 1, ids); +} + +template +void ImExPCTimestepper::buildOperators( + const std::map& ops, + const MHDFloat dt, const int size) +{ + std::map::const_iterator iOpA = + ops.find(ModelOperator::ImplicitLinear::id()); + std::map::const_iterator iOpB = + ops.find(ModelOperator::Time::id()); + std::map::const_iterator iOpC = + ops.find(ModelOperator::Boundary::id()); + + // Check if equation is solved in split form + bool isSplit = (ops.count(ModelOperator::SplitImplicitLinear::id()) > 0); + std::map::const_iterator iOpSC = + ops.find(ModelOperator::SplitBoundary::id()); + std::map::const_iterator iOpSA = + ops.find(ModelOperator::SplitImplicitLinear::id()); + std::map::const_iterator iOpSCV = + ops.find(ModelOperator::SplitBoundaryValue::id()); + + // Update timestep + this->mDt = dt; + + // Set explicit matrix + this->linearOperator(Tag::Operator::Rhs::id(), 0).resize(size, size); + details::addOperators(this->linearOperator(Tag::Operator::Rhs::id(), 0), 1.0, iOpA->second); + + // Set mass matrix + this->mMassMatrix.resize(size, size); + details::addOperators(this->mMassMatrix, 1.0, iOpB->second); + + // Set implicit matrix + for (int i = 0; i < this->steps(); ++i) + { + MHDFloat a = this->aIm(i); + + // Set LHS matrix + auto&& lhsMat = this->linearOperator(Tag::Operator::Lhs::id(), a); + lhsMat.resize(size, size); + details::addOperators(lhsMat, 1.0, iOpB->second); + details::addOperators(lhsMat, -a * this->mDt, iOpA->second); + details::addOperators(lhsMat, 1.0, iOpC->second); + + if (isSplit) + { + // Initialise influence matrices + MHDFloat aInf = 0.0; + this->initMatrices(Tag::Operator::Influence::id(), aInf); + + // Set other LHS matrix + auto&& infMatrix = + this->linearOperator(Tag::Operator::Influence::id(), aInf); + infMatrix.resize(size, size); + details::addOperators(infMatrix, 1.0, iOpSA->second); + details::addOperators(infMatrix, 1.0, iOpSC->second); + } + } + + if (isSplit) + { + // Store information for particular solution + auto&& infRhs = this->reg(Register::Influence::id()); + details::initInfluence(infRhs, iOpSCV->second, iOpSC->second); + } +} + +template +bool ImExPCTimestepper::preSolve() +{ + const bool hasInfluence = + this->hasLinearOperator(Tag::Operator::Influence::id()); + + // Use influence matrix + if (hasInfluence) + { + const bool isFirstPass = (this->mOpId == Tag::Operator::Lhs::id()); + if (isFirstPass) + { + this->mOpId = Tag::Operator::Influence::id(); + this->mId = 0.0; + + return true; + } + else + { + this->mOpId = Tag::Operator::Lhs::id(); + } + } + + // First step + if (this->mStep == 0) + { + // Reset error + if (this->mspScheme->useEmbedded()) + { + this->mError = 0.0; + } + + // Build RHS + MHDFloat aIm = (1.0 - this->mspScheme->aIm(this->mStep)) * this->mDt; + MHDFloat aMass = 1.0; + MHDFloat aN = this->mDt; + if (this->mHasExplicit) + { + details::computeSet(this->reg(Register::Explicit::id()), -1.0, + this->reg(Register::Rhs::id())); + details::computeSet(this->reg(Register::Rhs::id()), aN, + this->reg(Register::Explicit::id())); + } + details::computeAMXPY(this->reg(Register::Rhs::id()), + this->mMassMatrix, aMass, + this->reg(Register::Intermediate::id())); + details::computeAMXPY(this->reg(Register::Rhs::id()), + this->linearOperator(Tag::Operator::Rhs::id(), 0), aIm, + this->reg(Register::Intermediate::id())); + + this->mId = this->mspScheme->aIm(this->mStep); + + // Include inhomogeneous boundary conditions + this->addInhomogeneous(); + } + else + { + MHDFloat aNold = -this->mspScheme->aIm(this->mStep) * this->mDt; + MHDFloat aNnew = this->mspScheme->aIm(this->mStep) * this->mDt; + if (this->mHasExplicit) + { + details::computeSet(this->reg(Register::Error::id()), aNold, + this->reg(Register::Explicit::id())); + details::computeSet(this->reg(Register::Explicit::id()), -1.0, + this->reg(Register::Rhs::id())); + details::computeSet(this->reg(Register::Rhs::id()), aNnew, + this->reg(Register::Explicit::id())); + } + details::computeXPAY(this->reg(Register::Rhs::id()), + this->reg(Register::Error::id()), 1.0); + + this->mId = this->mspScheme->aIm(this->mStep); + } + + return true; +} + +template +bool ImExPCTimestepper::postSolve() +{ + const bool hasInfluence = + this->hasLinearOperator(Tag::Operator::Influence::id()); + + if (hasInfluence) + { + const bool isFirstPass = (this->mOpId == Tag::Operator::Influence::id()); + + if (isFirstPass) + { + // Apply quasi-inverse + details::computeMV(this->reg(Register::Rhs::id()), + this->mMassMatrix, + this->reg(Register::Solution::id())); + + return true; + } + else + { + // Correct solution with green's function + details::computeInfluenceCorrection( + this->reg(Register::Solution::id()), + this->reg(Register::Influence::id())); + } + } + + if (this->mStep == 0) + { + // Store predictor solution + details::computeSet(this->reg(Register::Intermediate::id()), + this->reg(Register::Solution::id())); + + this->mStep += 1; + } + else + { + // Build corrected solution + details::computeSet(this->reg(Register::Error::id()), + this->reg(Register::Solution::id())); + details::computeXPAY(this->reg(Register::Solution::id()), + this->reg(Register::Intermediate::id()), 1.0); + details::computeSet(this->reg(Register::Intermediate::id()), + this->reg(Register::Solution::id())); + + // Compute error + if (this->mspScheme->useEmbedded()) + { + details::computeErrorFromDiff(this->mError, + this->reg(Register::Error::id()), + this->reg(Register::Intermediate::id())); + } + + this->mStep += 1; + + // Check if we are done + if (this->mStep == this->steps()) + { + this->mStep = 0; + } + } + + return false; +} + +template +void ImExPCTimestepper::setSolution(const View::View>>& sol, const std::size_t startRow) +{ + details::computeSet(this->reg(Register::Solution::id()), + sol, startRow); +} + +template +void ImExPCTimestepper::addRhs(const View::View>>& rhs, const std::size_t startRow) +{ + details::computeAXPY(this->reg(Register::Rhs::id()), 1.0, + rhs, startRow); +} + +template +void ImExPCTimestepper::getSolution(View::View>>& sol, const std::size_t startRow) +{ + details::computeSet(sol, this->reg(Register::Solution::id()), startRow); +} + +} // namespace Views +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace QuICC + +#endif // QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_IMEXPCTIMSTEPPER_HPP diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/Tags.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/Tags.hpp new file mode 100644 index 00000000..0b3c7bdc --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/Tags.hpp @@ -0,0 +1,27 @@ +/** + * @file Tags.hpp + * @brief Tag types + */ +#pragma once + +// System includes +// + +// Project includes +// + +namespace QuICC { +namespace Timestep { +namespace PredictorCorrector { + +// +// Tags +// + +struct base_t +{ +}; + +} // namespace PredictorCorrector +} // namespace Timestep +} // namespace QuICC diff --git a/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/TimestepperCoordinator.hpp b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/TimestepperCoordinator.hpp new file mode 100644 index 00000000..db43779e --- /dev/null +++ b/Components/TimestepPcViews/Timestep/PredictorCorrector/Views/TimestepperCoordinator.hpp @@ -0,0 +1,503 @@ +/** + * @file SparseCoordinatorData.hpp + * @brief Implementation of the base for a general sparse solver coordinator + */ + +#ifndef QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_SPARSECOORDINATORDATA_HPP +#define QUICC_TIMESTEP_PREDICTORCORRECTOR_VIEWS_SPARSECOORDINATORDATA_HPP + +// System includes +// +#include + +// Project includes +// +#include "QuICC/Debug/DebuggerMacro.h" +#ifdef QUICC_DEBUG +#include "QuICC/PhysicalNames/Coordinator.hpp" +#include "QuICC/ModelOperator/Coordinator.hpp" +#include "QuICC/Tools/IdToHuman.hpp" +#endif +#include "Types/Typedefs.hpp" +#include "Timestep/PredictorCorrector/Views/Tags.hpp" +#include "View/ViewDense.hpp" + +#include +namespace QuICC { + +namespace Timestep { + +namespace PredictorCorrector { + +namespace Views { + + struct TimestepperInfo + { + bool isComplex; + std::size_t solverIndex; + std::size_t fieldIndex; + std::size_t rows; + std::size_t cols; + std::size_t blockN; + std::map ops; + }; + + template