From 437f57ecc297aaa194cfce3e0eb841cdf5a75a23 Mon Sep 17 00:00:00 2001 From: "Franz R. Sattler" Date: Wed, 18 Dec 2024 14:38:42 +0100 Subject: [PATCH 1/2] created a finite volume discretization --- .../discretization/FV/discretization.hh | 56 +++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 DiFfRG/include/DiFfRG/discretization/FV/discretization.hh diff --git a/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh b/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh new file mode 100644 index 0000000..741f39a --- /dev/null +++ b/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh @@ -0,0 +1,56 @@ +#pragma once + +// external libraries +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// DiFfRG +#include + +namespace DiFfRG +{ + namespace FV + { + using namespace dealii; + + /** + * @brief Class to manage the system on which we solve, i.e. fe spaces, grids, etc. + * This class is a System for CG systems. + * + * @tparam Model_ The Model class used for the Simulation + */ + template class Discretization + { + public: + using Components = Components_; + using NumberType = NumberType_; + using VectorType = Vector; + using SparseMatrixType = SparseMatrix; + using Mesh = Mesh_; + static constexpr uint dim = Mesh::dim; + + Discretization(Mesh &mesh, const JSONValue &json) : mesh(mesh), json(json) {}; + + const auto &get_mapping() const { return mapping; } + const auto &get_triangulation() const { return mesh.get_triangulation(); } + auto &get_triangulation() { return mesh.get_triangulation(); } + const auto &get_json() const { return json; } + + void reinit() {} + + protected: + Mesh &mesh; + JSONValue json; + + AffineConstraints constraints; + MappingQ1 mapping; + }; + } // namespace FV +} // namespace DiFfRG From 5610fd1bb06ea12c46599b1523c51cc4bd3e243a Mon Sep 17 00:00:00 2001 From: "Franz R. Sattler" Date: Sat, 21 Dec 2024 22:39:28 +0100 Subject: [PATCH 2/2] fix: bug in deal.ii, fixed by disabling taskflow in deal.ii fix: SolutionTransfer has now a different interface, adapted it feat: FlowingVariables for different discretizations !FlowingVariables->FE::FlowingVariables Signed-off-by: Franz R. Sattler --- CHANGELOG.md | 12 ++ DiFfRG/documentation/tutorials/tut1.md | 2 +- .../include/DiFfRG/discretization/FEM/cg.hh | 7 + .../include/DiFfRG/discretization/FEM/dg.hh | 7 + .../include/DiFfRG/discretization/FEM/ldg.hh | 7 + .../DiFfRG/discretization/data/data.hh | 197 ++++++++++++++---- .../discretization/mesh/h_adaptivity.hh | 2 +- .../DiFfRG/model/component_descriptor.hh | 7 - Examples/ONfiniteT/CG.cc | 2 +- Examples/ONfiniteT/LDG.cc | 2 +- Examples/ONfiniteT/dDG.cc | 2 +- Examples/ONfiniteT/parameter.json | 2 +- Examples/QuarkMesonLPAprime/.clang-format | 2 +- Examples/QuarkMesonLPAprime/QM.cc | 2 +- Examples/YangMills/Full/.clang-format | 2 +- Examples/YangMills/Full/Yang-Mills.cc | 3 +- Examples/YangMills/Full/model.hh | 93 +++++---- Examples/YangMills/SP/.clang-format | 2 +- Examples/YangMills/SP/Yang-Mills.cc | 3 +- Examples/YangMills/SP/model.hh | 52 +++-- README.md | 4 +- Tutorials/tut1/tut1.cc | 2 +- Tutorials/tut2/tut2.cc | 2 +- VERSION | 2 +- external/build_dealii.sh | 3 +- external/build_sundials.sh | 2 +- 26 files changed, 280 insertions(+), 143 deletions(-) create mode 100644 CHANGELOG.md diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..3ad2a0c --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,12 @@ +# Changelog + +## Version 1.0.1 + +### Fixed + +- There was a bug with TaskFlow internally in deal.ii. Fixed for now by simply disabling taskflow in deal.ii. +- deal.ii changed its interface for dealii::SolutionTransfer. Adapted the corresponding methods. + +### Changed + +- The FlowingVariables classes are now in separate namespaces. For finite elements, use DiFfRG::FE::FlowingVariables, for pure variable systems use DiFfRG::FlowingVariables. diff --git a/DiFfRG/documentation/tutorials/tut1.md b/DiFfRG/documentation/tutorials/tut1.md index ebfac8e..e51acb9 100644 --- a/DiFfRG/documentation/tutorials/tut1.md +++ b/DiFfRG/documentation/tutorials/tut1.md @@ -137,7 +137,7 @@ We now use the types we defined above to construct objects of all the classes de With everything prepared, we are now ready to set up and run the equation system. To do so, we create an initial condition on the discretization (i.e. FEM space) we set up before, ```Cpp // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); ``` and use it to run the time-stepper from RG-time 0 to the final RG-time, which we infer from the parameter file: diff --git a/DiFfRG/include/DiFfRG/discretization/FEM/cg.hh b/DiFfRG/include/DiFfRG/discretization/FEM/cg.hh index d4113c4..93d0d02 100644 --- a/DiFfRG/include/DiFfRG/discretization/FEM/cg.hh +++ b/DiFfRG/include/DiFfRG/discretization/FEM/cg.hh @@ -80,6 +80,13 @@ namespace DiFfRG void reinit() { setup_dofs(); } + std::vector get_block_structure() const + { + std::vector block_structure{dof_handler.n_dofs()}; + if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables()); + return block_structure; + } + protected: void setup_dofs() { diff --git a/DiFfRG/include/DiFfRG/discretization/FEM/dg.hh b/DiFfRG/include/DiFfRG/discretization/FEM/dg.hh index 2303968..951d093 100644 --- a/DiFfRG/include/DiFfRG/discretization/FEM/dg.hh +++ b/DiFfRG/include/DiFfRG/discretization/FEM/dg.hh @@ -94,6 +94,13 @@ namespace DiFfRG return dof; } + std::vector get_block_structure() const + { + std::vector block_structure{dof_handler.n_dofs()}; + if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables()); + return block_structure; + } + protected: void setup_dofs() { diff --git a/DiFfRG/include/DiFfRG/discretization/FEM/ldg.hh b/DiFfRG/include/DiFfRG/discretization/FEM/ldg.hh index 662fbc0..6b16628 100644 --- a/DiFfRG/include/DiFfRG/discretization/FEM/ldg.hh +++ b/DiFfRG/include/DiFfRG/discretization/FEM/ldg.hh @@ -88,6 +88,13 @@ namespace DiFfRG return dof; } + std::vector get_block_structure() const + { + std::vector block_structure{dof_handler[0]->n_dofs()}; + if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables()); + return block_structure; + } + protected: void setup_dofs() { diff --git a/DiFfRG/include/DiFfRG/discretization/data/data.hh b/DiFfRG/include/DiFfRG/discretization/data/data.hh index 89d36b5..108dee2 100644 --- a/DiFfRG/include/DiFfRG/discretization/data/data.hh +++ b/DiFfRG/include/DiFfRG/discretization/data/data.hh @@ -27,21 +27,6 @@ namespace DiFfRG private: FUN fun; }; - - struct NoDiscretization { - static constexpr uint dim = 0; - using NumberType = double; - using Components = void; - - const DoFHandler &get_dof_handler() const { throw std::runtime_error("No discretization chosen!"); } - - template static std::vector get_block_structure(const Model &) - { - std::vector block_structure{0}; - block_structure.push_back(Model::Components::count_variables()); - return block_structure; - } - }; } // namespace internal /** @@ -51,29 +36,15 @@ namespace DiFfRG * * @tparam Discretization Spatial Discretization used in the system */ - template - class FlowingVariables : public AbstractFlowingVariables + template class FlowingVariables : public AbstractFlowingVariables { public: - using NumberType = typename Discretization::NumberType; - using Components = typename Discretization::Components; - static constexpr uint dim = Discretization::dim; + using NumberType = NT; /** * @brief Construct a new Flowing Variables object - * - * @param discretization The spatial discretization to use */ - FlowingVariables(const Discretization &discretization) : dof_handler(discretization.get_dof_handler()) - { - auto block_structure = Components::get_block_structure(dof_handler); - m_data = (block_structure); - } - - /** - * @brief Construct a new Flowing Variables object - */ - FlowingVariables() : dof_handler(*static_cast *>(0)) {} + FlowingVariables() {} /** * @brief Interpolates the initial condition from a numerical model. @@ -83,15 +54,9 @@ namespace DiFfRG */ template void interpolate(const Model &model) { - if constexpr (Model::Components::count_fe_functions() > 0) { - auto interpolating_function = [&model](const auto &p, auto &values) { model.initial_condition(p, values); }; - internal::FunctionFromLambda initial_condition_function(std::move(interpolating_function), - dof_handler.get_fe().n_components()); - VectorTools::interpolate(dof_handler, initial_condition_function, m_data.block(0)); - } else { - auto block_structure = Discretization::get_block_structure(model); - m_data = (block_structure); - } + std::vector block_structure{0}; + block_structure.push_back(Model::Components::count_variables()); + m_data = (block_structure); if (m_data.n_blocks() > 1) model.initial_condition_variables(m_data.block(1)); } @@ -120,8 +85,156 @@ namespace DiFfRG virtual const Vector &variable_data() const override { return m_data.block(1); } private: - const DoFHandler &dof_handler; BlockVector m_data; }; + namespace FE + { + /** + * @brief A class to set up initial data for whatever discretization we have chosen. + * Also used to switch/manage memory, vectors, matrices over interfaces between spatial discretization and + * separate variables. + * + * @tparam Discretization Spatial Discretization used in the system + */ + template + class FlowingVariables : public AbstractFlowingVariables + { + public: + using NumberType = typename Discretization::NumberType; + using Components = typename Discretization::Components; + static constexpr uint dim = Discretization::dim; + + /** + * @brief Construct a new Flowing Variables object + * + * @param discretization The spatial discretization to use + */ + FlowingVariables(const Discretization &discretization) + : discretization(discretization), dof_handler(discretization.get_dof_handler()) + { + } + + /** + * @brief Interpolates the initial condition from a numerical model. + * + * @param model The model to interpolate from. Must provide a method initial_condition(const Point &, + * Vector &) + */ + template void interpolate(const Model &model) + { + auto block_structure = discretization.get_block_structure(); + m_data = (block_structure); + + if constexpr (Model::Components::count_fe_functions() > 0) { + auto interpolating_function = [&model](const auto &p, auto &values) { model.initial_condition(p, values); }; + internal::FunctionFromLambda initial_condition_function(std::move(interpolating_function), + dof_handler.get_fe().n_components()); + VectorTools::interpolate(dof_handler, initial_condition_function, m_data.block(0)); + } + if (m_data.n_blocks() > 1) model.initial_condition_variables(m_data.block(1)); + } + + /** + * @brief Obtain the data vector holding both spatial (block 0) and variable (block 1) data. + * + * @return BlockVector& The data vector. + */ + virtual BlockVector &data() override { return m_data; } + virtual const BlockVector &data() const override { return m_data; } + + /** + * @brief Obtain the spatial data vector. + * + * @return Vector& The spatial data vector. + */ + virtual Vector &spatial_data() override { return m_data.block(0); } + virtual const Vector &spatial_data() const override { return m_data.block(0); } + + /** + * @brief Obtain the variable data vector. + * + * @return Vector& The variable data vector. + */ + virtual Vector &variable_data() override { return m_data.block(1); } + virtual const Vector &variable_data() const override { return m_data.block(1); } + + private: + const Discretization &discretization; + const DoFHandler &dof_handler; + BlockVector m_data; + }; + } // namespace FE + + namespace FV + { + /** + * @brief A class to set up initial data for whatever discretization we have chosen. + * Also used to switch/manage memory, vectors, matrices over interfaces between spatial discretization and + * separate variables. + * + * @tparam Discretization Spatial Discretization used in the system + */ + template + class FlowingVariables : public AbstractFlowingVariables + { + public: + using NumberType = typename Discretization::NumberType; + using Components = typename Discretization::Components; + static constexpr uint dim = Discretization::dim; + + /** + * @brief Construct a new Flowing Variables object + * + * @param discretization The spatial discretization to use + */ + FlowingVariables(const Discretization &discretization) : discretization(discretization) {} + + /** + * @brief Interpolates the initial condition from a numerical model. + * + * @param model The model to interpolate from. Must provide a method initial_condition(const Point &, + * Vector &) + */ + template void interpolate(const Model &model) + { + auto block_structure = discretization.get_block_structure(); + m_data = (block_structure); + + if constexpr (Model::Components::count_fe_functions() > 0) { + // TODO + } + if (m_data.n_blocks() > 1) model.initial_condition_variables(m_data.block(1)); + } + + /** + * @brief Obtain the data vector holding both spatial (block 0) and variable (block 1) data. + * + * @return BlockVector& The data vector. + */ + virtual BlockVector &data() override { return m_data; } + virtual const BlockVector &data() const override { return m_data; } + + /** + * @brief Obtain the spatial data vector. + * + * @return Vector& The spatial data vector. + */ + virtual Vector &spatial_data() override { return m_data.block(0); } + virtual const Vector &spatial_data() const override { return m_data.block(0); } + + /** + * @brief Obtain the variable data vector. + * + * @return Vector& The variable data vector. + */ + virtual Vector &variable_data() override { return m_data.block(1); } + virtual const Vector &variable_data() const override { return m_data.block(1); } + + private: + const Discretization &discretization; + BlockVector m_data; + }; + } // namespace FV + } // namespace DiFfRG diff --git a/DiFfRG/include/DiFfRG/discretization/mesh/h_adaptivity.hh b/DiFfRG/include/DiFfRG/discretization/mesh/h_adaptivity.hh index a8be596..b62785b 100644 --- a/DiFfRG/include/DiFfRG/discretization/mesh/h_adaptivity.hh +++ b/DiFfRG/include/DiFfRG/discretization/mesh/h_adaptivity.hh @@ -101,7 +101,7 @@ namespace DiFfRG assembler.reinit(); assembler.reinit_vector(solution); - solution_trans.interpolate(previous_solution, solution); + solution_trans.interpolate(solution); constraints.distribute(solution); return true; diff --git a/DiFfRG/include/DiFfRG/model/component_descriptor.hh b/DiFfRG/include/DiFfRG/model/component_descriptor.hh index 84f11db..a88ffde 100644 --- a/DiFfRG/include/DiFfRG/model/component_descriptor.hh +++ b/DiFfRG/include/DiFfRG/model/component_descriptor.hh @@ -165,13 +165,6 @@ namespace DiFfRG return dofs_per_component; } - template static std::vector get_block_structure(const DoFH &dofh) - { - std::vector block_structure{dofh.n_dofs()}; - if (n_variables > 0) block_structure.push_back(n_variables); - return block_structure; - } - private: template using SubsystemMatrix = std::array, n_fe_subsystems>; diff --git a/Examples/ONfiniteT/CG.cc b/Examples/ONfiniteT/CG.cc index 75074db..f2232f9 100644 --- a/Examples/ONfiniteT/CG.cc +++ b/Examples/ONfiniteT/CG.cc @@ -33,7 +33,7 @@ int main(int argc, char *argv[]) TimeStepper time_stepper(json, &assembler, &data_out, &mesh_adaptor); // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); // Now we start the timestepping diff --git a/Examples/ONfiniteT/LDG.cc b/Examples/ONfiniteT/LDG.cc index 8e269eb..8091d6e 100644 --- a/Examples/ONfiniteT/LDG.cc +++ b/Examples/ONfiniteT/LDG.cc @@ -33,7 +33,7 @@ int main(int argc, char *argv[]) TimeStepper time_stepper(json, &assembler, &data_out, &mesh_adaptor); // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); // Now we start the timestepping diff --git a/Examples/ONfiniteT/dDG.cc b/Examples/ONfiniteT/dDG.cc index 7ba030b..2e3aa7f 100644 --- a/Examples/ONfiniteT/dDG.cc +++ b/Examples/ONfiniteT/dDG.cc @@ -33,7 +33,7 @@ int main(int argc, char *argv[]) TimeStepper time_stepper(json, &assembler, &data_out, &mesh_adaptor); // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); // Now we start the timestepping diff --git a/Examples/ONfiniteT/parameter.json b/Examples/ONfiniteT/parameter.json index 472dbe8..29cb7d4 100644 --- a/Examples/ONfiniteT/parameter.json +++ b/Examples/ONfiniteT/parameter.json @@ -60,7 +60,7 @@ } }, "output": { - "verbosity": 0, + "verbosity": 1, "folder": "./", "name": "output" } diff --git a/Examples/QuarkMesonLPAprime/.clang-format b/Examples/QuarkMesonLPAprime/.clang-format index ca4d744..75085d1 100644 --- a/Examples/QuarkMesonLPAprime/.clang-format +++ b/Examples/QuarkMesonLPAprime/.clang-format @@ -5,7 +5,7 @@ TabWidth: 2 BreakBeforeBraces: Linux AllowShortIfStatementsOnASingleLine: true IndentCaseLabels: false -ColumnLimit: 180 +ColumnLimit: 120 AccessModifierOffset: -2 NamespaceIndentation: All AllowShortEnumsOnASingleLine: true diff --git a/Examples/QuarkMesonLPAprime/QM.cc b/Examples/QuarkMesonLPAprime/QM.cc index 08c0f41..87fdca2 100644 --- a/Examples/QuarkMesonLPAprime/QM.cc +++ b/Examples/QuarkMesonLPAprime/QM.cc @@ -32,7 +32,7 @@ int main(int argc, char *argv[]) TimeStepper time_stepper(json, &assembler, &data_out); // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); // Now we start the timestepping diff --git a/Examples/YangMills/Full/.clang-format b/Examples/YangMills/Full/.clang-format index ca4d744..75085d1 100644 --- a/Examples/YangMills/Full/.clang-format +++ b/Examples/YangMills/Full/.clang-format @@ -5,7 +5,7 @@ TabWidth: 2 BreakBeforeBraces: Linux AllowShortIfStatementsOnASingleLine: true IndentCaseLabels: false -ColumnLimit: 180 +ColumnLimit: 120 AccessModifierOffset: -2 NamespaceIndentation: All AllowShortEnumsOnASingleLine: true diff --git a/Examples/YangMills/Full/Yang-Mills.cc b/Examples/YangMills/Full/Yang-Mills.cc index a517dcf..3edf47d 100644 --- a/Examples/YangMills/Full/Yang-Mills.cc +++ b/Examples/YangMills/Full/Yang-Mills.cc @@ -25,8 +25,7 @@ bool run(const JSONValue &json, const std::string logger) // Define the objects needed to run the simulation Model model(json); Assembler assembler(model, json); - DataOutput<0, VectorType> data_out(config_helper.get_top_folder(), config_helper.get_output_name(), config_helper.get_output_folder(), json); - TimeStepper time_stepper(&assembler, data_out, nullptr, json); + TimeStepper time_stepper(json, &assembler); // Set up the initial condition FlowingVariables initial_condition; diff --git a/Examples/YangMills/Full/model.hh b/Examples/YangMills/Full/model.hh index baf46aa..51e1fb6 100644 --- a/Examples/YangMills/Full/model.hh +++ b/Examples/YangMills/Full/model.hh @@ -21,7 +21,7 @@ struct Parameters { tilt_A3 = json.get_double("/physical/tilt_A3"); tilt_A4 = json.get_double("/physical/tilt_A4"); tilt_Acbc = json.get_double("/physical/tilt_Acbc"); - + m2A = json.get_double("/physical/m2A"); p_grid_min = json.get_double("/discretization/p_grid_min"); @@ -45,23 +45,21 @@ static constexpr uint p_grid_size = 96; static constexpr uint S1_size = 8; static constexpr uint SPhi_size = 7; // As for components, we have one FE function (u) and several extractors. -using VariableDesc = VariableDescriptor< - Scalar<"m2A">, +using VariableDesc = + VariableDescriptor, - FunctionND<"ZA3", p_grid_size, S1_size, SPhi_size>, - FunctionND<"ZAcbc", p_grid_size, S1_size, SPhi_size>, - FunctionND<"ZA4tadpole", p_grid_size, S1_size, SPhi_size>, - FunctionND<"ZA4SP", p_grid_size>, + FunctionND<"ZA3", p_grid_size, S1_size, SPhi_size>, + FunctionND<"ZAcbc", p_grid_size, S1_size, SPhi_size>, + FunctionND<"ZA4tadpole", p_grid_size, S1_size, SPhi_size>, FunctionND<"ZA4SP", p_grid_size>, - FunctionND<"ZA", p_grid_size>, - FunctionND<"Zc", p_grid_size> - >; + FunctionND<"ZA", p_grid_size>, FunctionND<"Zc", p_grid_size>>; using Components = ComponentDescriptor, VariableDesc, ExtractorDescriptor<>>; constexpr auto idxv = VariableDesc{}; /** - * @brief This class implements the numerical model for the quark-meson model at finite temperature and chemical potential. + * @brief This class implements the numerical model for the quark-meson model at finite temperature and chemical + * potential. */ class YangMills : public def::AbstractModel, public def::fRG, // this handles the fRG time @@ -71,7 +69,8 @@ class YangMills : public def::AbstractModel, using Coordinates1D = LogarithmicCoordinates1D; using CoordinatesAng = LinearCoordinates1D; - using Coordinates3D = CoordinatePackND, LinearCoordinates1D, LinearCoordinates1D>; + using Coordinates3D = + CoordinatePackND, LinearCoordinates1D, LinearCoordinates1D>; const Coordinates1D coordinates1D; const CoordinatesAng S1_coordinates; const CoordinatesAng SPhi_coordinates; @@ -89,18 +88,14 @@ class YangMills : public def::AbstractModel, public: YangMills(const JSONValue &json) - : def::fRG(json.get_double("/physical/Lambda")), prm(json), - coordinates1D(p_grid_size, prm.p_grid_min, prm.p_grid_max, prm.p_grid_bias), - S1_coordinates(S1_size, 0.0, 0.9999), - SPhi_coordinates(SPhi_size, 0.0, 2. * M_PI), - coordinates3D(coordinates1D, S1_coordinates, SPhi_coordinates), - grid1D(make_grid(coordinates1D)), - S1_grid(make_grid(S1_coordinates)), - SPhi_grid(make_grid(SPhi_coordinates)), - grid3D(make_grid(coordinates3D)), - flow_equations(json), - dtZc(coordinates1D), dtZA(coordinates1D), ZA(coordinates1D), Zc(coordinates1D), // propagators - ZA4SP(coordinates1D), // 1D couplings + : def::fRG(json.get_double("/physical/Lambda")), prm(json), + coordinates1D(p_grid_size, prm.p_grid_min, prm.p_grid_max, prm.p_grid_bias), + S1_coordinates(S1_size, 0.0, 0.9999), SPhi_coordinates(SPhi_size, 0.0, 2. * M_PI), + coordinates3D(coordinates1D, S1_coordinates, SPhi_coordinates), grid1D(make_grid(coordinates1D)), + S1_grid(make_grid(S1_coordinates)), SPhi_grid(make_grid(SPhi_coordinates)), grid3D(make_grid(coordinates3D)), + flow_equations(json), dtZc(coordinates1D), dtZA(coordinates1D), ZA(coordinates1D), + Zc(coordinates1D), // propagators + ZA4SP(coordinates1D), // 1D couplings ZAcbc(coordinates3D), ZA3(coordinates3D), ZA4tadpole(coordinates3D) // 3D couplings { flow_equations.set_k(prm.Lambda); @@ -112,9 +107,12 @@ public: for (uint i = 0; i < p_grid_size; ++i) { for (uint j = 0; j < S1_size; ++j) { for (uint k = 0; k < SPhi_size; ++k) { - values[idxv("ZA3") + i * S1_size * SPhi_size + j * SPhi_size + k] = std::sqrt(4. * M_PI * prm.alphaA3) + prm.tilt_A3 * std::log(grid1D[i] / prm.p_grid_max); - values[idxv("ZAcbc") + i * S1_size * SPhi_size + j * SPhi_size + k] = std::sqrt(4. * M_PI * prm.alphaAcbc) + prm.tilt_Acbc * std::log(grid1D[i] / prm.p_grid_max); - values[idxv("ZA4tadpole") + i * S1_size * SPhi_size + j * SPhi_size + k] = 4. * M_PI * prm.alphaA4 + prm.tilt_A4 * std::log(grid1D[i] / prm.p_grid_max); + values[idxv("ZA3") + i * S1_size * SPhi_size + j * SPhi_size + k] = + std::sqrt(4. * M_PI * prm.alphaA3) + prm.tilt_A3 * std::log(grid1D[i] / prm.p_grid_max); + values[idxv("ZAcbc") + i * S1_size * SPhi_size + j * SPhi_size + k] = + std::sqrt(4. * M_PI * prm.alphaAcbc) + prm.tilt_Acbc * std::log(grid1D[i] / prm.p_grid_max); + values[idxv("ZA4tadpole") + i * S1_size * SPhi_size + j * SPhi_size + k] = + 4. * M_PI * prm.alphaA4 + prm.tilt_A4 * std::log(grid1D[i] / prm.p_grid_max); } } @@ -147,11 +145,9 @@ public: ZA.update(&variables.data()[idxv("ZA")]); Zc.update(&variables.data()[idxv("Zc")]); - // set up arguments for the integrators - const auto arguments = std::tie(ZA3, ZAcbc, ZA4SP, ZA4tadpole, - dtZc, Zc, dtZA, ZA, m2A); + const auto arguments = std::tie(ZA3, ZAcbc, ZA4SP, ZA4tadpole, dtZc, Zc, dtZA, ZA, m2A); // copy the propagators for comparison std::vector old_dtZA(p_grid_size); @@ -207,9 +203,11 @@ public: residual[idxv("Zc") + i] = dtZc.data()[i]; } - std::apply([&](const auto &...args) { residual[idxv("m2A")] = flow_equations.m2A_integrator.get(k, 0., args...); }, arguments); + std::apply( + [&](const auto &...args) { residual[idxv("m2A")] = flow_equations.m2A_integrator.get(k, 0., args...); }, + arguments); - if(Zc[0] < 0) throw std::runtime_error("Zc < 0"); + if (Zc[0] < 0) throw std::runtime_error("Zc < 0"); /*for (uint i = 0; i < p_grid_size; ++i) { if(!std::isfinite(residual[idxv("ZA") + i])) std::cout << "res ZA not finite" << std::endl; @@ -217,8 +215,9 @@ public: for (uint j = 0; j < S1_size; ++j) { for (uint k = 0; k < SPhi_size; ++k) { - if(!std::isfinite(residual[idxv("ZA3") + i * S1_size * SPhi_size + j * SPhi_size + k])) std::cout << "res ZA3 not finite" << std::endl; - if(!std::isfinite(residual[idxv("ZAcbc") + i * S1_size * SPhi_size + j * SPhi_size + k])) std::cout << "res ZAcbc not finite" << std::endl; + if(!std::isfinite(residual[idxv("ZA3") + i * S1_size * SPhi_size + j * SPhi_size + k])) std::cout << "res ZA3 + not finite" << std::endl; if(!std::isfinite(residual[idxv("ZAcbc") + i * S1_size * SPhi_size + j * SPhi_size + k])) + std::cout << "res ZAcbc not finite" << std::endl; } } @@ -226,7 +225,8 @@ public: }*/ } - template void readouts(DataOut &output, const Point &, const Solutions &sol) const + template + void readouts(DataOut &output, const Point &, const Solutions &sol) const { const auto &variables = get<"variables">(sol); const auto m2A = variables[idxv("m2A")]; @@ -244,7 +244,7 @@ public: Zs_data[i][0] = k; Zs_data[i][1] = grid1D[i]; Zs_data[i][2] = ZA[i]; - Zs_data[i][3] = (ZA[i]);// * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); + Zs_data[i][3] = (ZA[i]); // * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); Zs_data[i][4] = Zc[i]; } @@ -253,7 +253,7 @@ public: const auto *ZA4SP = &variables.data()[idxv("ZA4SP")]; std::vector> strong_couplings_data(grid1D.size(), std::vector(5, 0.)); for (uint i = 0; i < p_grid_size; ++i) { - const double ZA_p = (ZA[i]);// * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); + const double ZA_p = (ZA[i]); // * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); const double Zc_p = Zc[i]; const size_t cur_idx = i * S1_size * SPhi_size; strong_couplings_data[i][0] = k; @@ -274,16 +274,16 @@ public: const auto S1 = S1_grid[j]; const auto SPhi = SPhi_grid[k]; - const auto p = S0*(sqrt(1.f-1.f*S1*sin(SPhi))); - const auto r = S0*(sqrt(1.f+0.5f*S1*(1.7320508075688772f*cos(SPhi)+sin(SPhi)))); - const auto s = 0.7071067811865475f*S0*(sqrt(2.f-1.7320508075688772f*S1*cos(SPhi)+S1*sin(SPhi))); + const auto p = S0 * (sqrt(1.f - 1.f * S1 * sin(SPhi))); + const auto r = S0 * (sqrt(1.f + 0.5f * S1 * (1.7320508075688772f * cos(SPhi) + sin(SPhi)))); + const auto s = 0.7071067811865475f * S0 * (sqrt(2.f - 1.7320508075688772f * S1 * cos(SPhi) + S1 * sin(SPhi))); const double Zc_r = this->Zc(r); const double Zc_s = this->Zc(s); - const double Z_p = (this->ZA(p));// * powr<2>(p) + m2A) / powr<2>(p); - const double Z_r = (this->ZA(r));// * powr<2>(r) + m2A) / powr<2>(r); - const double Z_s = (this->ZA(s));// * powr<2>(s) + m2A) / powr<2>(s); + const double Z_p = (this->ZA(p)); // * powr<2>(p) + m2A) / powr<2>(p); + const double Z_r = (this->ZA(r)); // * powr<2>(r) + m2A) / powr<2>(r); + const double Z_s = (this->ZA(s)); // * powr<2>(s) + m2A) / powr<2>(s); strong_couplings_3D_data[cur_idx][0] = this->k; strong_couplings_3D_data[cur_idx][1] = grid3D[cur_idx][0]; @@ -299,8 +299,10 @@ public: } if (is_close(t, 0.)) { - const std::vector strong_couplings_header = std::vector{"k [GeV]", "p [GeV]", "alphaAcbc", "alphaA3", "alphaA4"}; - const std::vector strong_couplings_3D_header = std::vector{"k [GeV]", "S0 [GeV]", "S1", "SPhi", "ZAcbc", "ZA3", "ZA4tadpole", "alphaA3", "alphaAcbc"}; + const std::vector strong_couplings_header = + std::vector{"k [GeV]", "p [GeV]", "alphaAcbc", "alphaA3", "alphaA4"}; + const std::vector strong_couplings_3D_header = std::vector{ + "k [GeV]", "S0 [GeV]", "S1", "SPhi", "ZAcbc", "ZA3", "ZA4tadpole", "alphaA3", "alphaAcbc"}; const std::vector Zs_header = std::vector{"k [GeV]", "p [GeV]", "ZAbar", "ZA", "Zc"}; output.dump_to_csv("strong_couplings.csv", strong_couplings_data, false, strong_couplings_header); @@ -311,6 +313,5 @@ public: output.dump_to_csv("strong_couplings_3D.csv", strong_couplings_3D_data, true); output.dump_to_csv("Zs.csv", Zs_data, true); } - } }; diff --git a/Examples/YangMills/SP/.clang-format b/Examples/YangMills/SP/.clang-format index ca4d744..75085d1 100644 --- a/Examples/YangMills/SP/.clang-format +++ b/Examples/YangMills/SP/.clang-format @@ -5,7 +5,7 @@ TabWidth: 2 BreakBeforeBraces: Linux AllowShortIfStatementsOnASingleLine: true IndentCaseLabels: false -ColumnLimit: 180 +ColumnLimit: 120 AccessModifierOffset: -2 NamespaceIndentation: All AllowShortEnumsOnASingleLine: true diff --git a/Examples/YangMills/SP/Yang-Mills.cc b/Examples/YangMills/SP/Yang-Mills.cc index 40f2c97..14b8027 100644 --- a/Examples/YangMills/SP/Yang-Mills.cc +++ b/Examples/YangMills/SP/Yang-Mills.cc @@ -25,8 +25,7 @@ bool run(const JSONValue &json, const std::string logger) // Define the objects needed to run the simulation Model model(json); Assembler assembler(model, json); - DataOutput<0, VectorType> data_out(config_helper.get_top_folder(), config_helper.get_output_name(), config_helper.get_output_folder(), json); - TimeStepper time_stepper(&assembler, data_out, nullptr, json); + TimeStepper time_stepper(json, &assembler); // Set up the initial condition FlowingVariables initial_condition; diff --git a/Examples/YangMills/SP/model.hh b/Examples/YangMills/SP/model.hh index 1383f53..3b8df33 100644 --- a/Examples/YangMills/SP/model.hh +++ b/Examples/YangMills/SP/model.hh @@ -21,7 +21,7 @@ struct Parameters { tilt_A3 = json.get_double("/physical/tilt_A3"); tilt_A4 = json.get_double("/physical/tilt_A4"); tilt_Acbc = json.get_double("/physical/tilt_Acbc"); - + m2A = json.get_double("/physical/m2A"); p_grid_min = json.get_double("/discretization/p_grid_min"); @@ -43,22 +43,19 @@ struct Parameters { // Size of the momentum grid static constexpr uint p_grid_size = 96; // As for components, we have one FE function (u) and several extractors. -using VariableDesc = VariableDescriptor< - Scalar<"m2A">, +using VariableDesc = + VariableDescriptor, - FunctionND<"ZA3", p_grid_size>, - FunctionND<"ZAcbc", p_grid_size>, - FunctionND<"ZA4", p_grid_size>, + FunctionND<"ZA3", p_grid_size>, FunctionND<"ZAcbc", p_grid_size>, FunctionND<"ZA4", p_grid_size>, - FunctionND<"ZA", p_grid_size>, - FunctionND<"Zc", p_grid_size> - >; + FunctionND<"ZA", p_grid_size>, FunctionND<"Zc", p_grid_size>>; using Components = ComponentDescriptor, VariableDesc, ExtractorDescriptor<>>; constexpr auto idxv = VariableDesc{}; /** - * @brief This class implements the numerical model for the quark-meson model at finite temperature and chemical potential. + * @brief This class implements the numerical model for the quark-meson model at finite temperature and chemical + * potential. */ class YangMills : public def::AbstractModel, public def::fRG, // this handles the fRG time @@ -74,16 +71,15 @@ class YangMills : public def::AbstractModel, mutable YangMillsFlowEquations flow_equations; mutable TexLinearInterpolator1D dtZc, dtZA, ZA, Zc; - mutable TexLinearInterpolator1D ZAcbc, ZA3, ZA4; + mutable TexLinearInterpolator1D ZA4, ZAcbc, ZA3; public: YangMills(const JSONValue &json) - : def::fRG(json.get_double("/physical/Lambda")), prm(json), - coordinates1D(p_grid_size, prm.p_grid_min, prm.p_grid_max, prm.p_grid_bias), - grid1D(make_grid(coordinates1D)), - flow_equations(json), - dtZc(coordinates1D), dtZA(coordinates1D), ZA(coordinates1D), Zc(coordinates1D), // propagators - ZA4(coordinates1D), ZAcbc(coordinates1D), ZA3(coordinates1D) // couplings + : def::fRG(json.get_double("/physical/Lambda")), prm(json), + coordinates1D(p_grid_size, prm.p_grid_min, prm.p_grid_max, prm.p_grid_bias), grid1D(make_grid(coordinates1D)), + flow_equations(json), dtZc(coordinates1D), dtZA(coordinates1D), ZA(coordinates1D), + Zc(coordinates1D), // propagators + ZA4(coordinates1D), ZAcbc(coordinates1D), ZA3(coordinates1D) // couplings { flow_equations.set_k(prm.Lambda); flow_equations.print_parameters("log"); @@ -94,7 +90,8 @@ public: for (uint i = 0; i < p_grid_size; ++i) { values[idxv("ZA4") + i] = 4. * M_PI * prm.alphaA4 + prm.tilt_A4 * std::log(grid1D[i] / prm.p_grid_max); values[idxv("ZA3") + i] = std::sqrt(4. * M_PI * prm.alphaA3) + prm.tilt_A3 * std::log(grid1D[i] / prm.p_grid_max); - values[idxv("ZAcbc") + i] = std::sqrt(4. * M_PI * prm.alphaAcbc) + prm.tilt_Acbc * std::log(grid1D[i] / prm.p_grid_max); + values[idxv("ZAcbc") + i] = + std::sqrt(4. * M_PI * prm.alphaAcbc) + prm.tilt_Acbc * std::log(grid1D[i] / prm.p_grid_max); values[idxv("ZA") + i] = (powr<2>(grid1D[i]) + prm.m2A) / powr<2>(grid1D[i]); values[idxv("Zc") + i] = 1.; @@ -122,11 +119,9 @@ public: ZA.update(&variables.data()[idxv("ZA")]); Zc.update(&variables.data()[idxv("Zc")]); - // set up arguments for the integrators - const auto arguments = std::tie(ZA3, ZAcbc, ZA4, - dtZc, Zc, dtZA, ZA, m2A); + const auto arguments = std::tie(ZA3, ZAcbc, ZA4, dtZc, Zc, dtZA, ZA, m2A); // copy the propagators for comparison std::vector old_dtZA(p_grid_size); @@ -180,10 +175,13 @@ public: residual[idxv("Zc") + i] = dtZc.data()[i]; } - std::apply([&](const auto &...args) { residual[idxv("m2A")] = flow_equations.m2A_integrator.get(k, 0., args...); }, arguments); + std::apply( + [&](const auto &...args) { residual[idxv("m2A")] = flow_equations.m2A_integrator.get(k, 0., args...); }, + arguments); } - template void readouts(DataOut &output, const Point &, const Solutions &sol) const + template + void readouts(DataOut &output, const Point &, const Solutions &sol) const { const auto &variables = get<"variables">(sol); const auto m2A = variables[idxv("m2A")]; @@ -201,7 +199,7 @@ public: Zs_data[i][0] = k; Zs_data[i][1] = grid1D[i]; Zs_data[i][2] = ZA[i]; - Zs_data[i][3] = (ZA[i]);// * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); + Zs_data[i][3] = (ZA[i]); // * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); Zs_data[i][4] = Zc[i]; } @@ -210,7 +208,7 @@ public: const auto *ZA4 = &variables.data()[idxv("ZA4")]; std::vector> strong_couplings_data(grid1D.size(), std::vector(8, 0.)); for (uint i = 0; i < p_grid_size; ++i) { - const double ZA_p = (ZA[i]);// * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); + const double ZA_p = (ZA[i]); // * powr<2>(grid1D[i]) + m2A) / powr<2>(grid1D[i]); const double Zc_p = Zc[i]; strong_couplings_data[i][0] = k; strong_couplings_data[i][1] = grid1D[i]; @@ -223,7 +221,8 @@ public: } if (is_close(t, 0.)) { - const std::vector strong_couplings_header = std::vector{"k [GeV]", "p [GeV]", "alphaAcbc", "alphaA3", "alphaA4", "ZAcbc", "ZA3", "ZA4"}; + const std::vector strong_couplings_header = + std::vector{"k [GeV]", "p [GeV]", "alphaAcbc", "alphaA3", "alphaA4", "ZAcbc", "ZA3", "ZA4"}; const std::vector Zs_header = std::vector{"k [GeV]", "p [GeV]", "ZAbar", "ZA", "Zc"}; output.dump_to_csv("strong_couplings.csv", strong_couplings_data, false, strong_couplings_header); @@ -232,6 +231,5 @@ public: output.dump_to_csv("strong_couplings.csv", strong_couplings_data, true); output.dump_to_csv("Zs.csv", Zs_data, true); } - } }; diff --git a/README.md b/README.md index 244fb18..5eb15d1 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ [![arXiv](https://img.shields.io/badge/arXiv-2412.13043-b31b1b.svg?style=for-the-badge)](https://arxiv.org/abs/2412.13043) -[![Doxygen](https://img.shields.io/badge/doxygen-2C4AA8?style=for-the-badge&logo=doxygen&logoColor=white)](https://satfra.github.io/DiFfRG/cpp/index.html) +[![Doxygen](https://img.shields.io/badge/doxygen-2C4AA8?style=for-the-badge&logo=c%2B%2B&logoColor=white)](https://satfra.github.io/DiFfRG/cpp/index.html) [![Wolfram](https://img.shields.io/badge/wolfram_doc-cf1c10?style=for-the-badge&logo=wolfram)](https://satfra.github.io/DiFfRG/wolfram/html/guide/DiFfRG.html) [![Python](https://img.shields.io/badge/python_doc-3670A0?style=for-the-badge&logo=python&logoColor=ffdd54)](https://satfra.github.io/DiFfRG/python/index.html) @@ -21,7 +21,7 @@ Both explicit and implicit timestepping methods are available and allow thus for We also include a set of tools for the evaluation of integrals and discretization of momentum dependencies. -For an overview, please see the [accompanying paper](https://arxiv.org/abs/...), the ***[tutorial page](https://satfra.github.io/DiFfRG/cpp/TutorialTOC.html)*** in the [documentation](https://satfra.github.io/DiFfRG/cpp/index.html) and the examples in `Examples/`. +For an overview, please see the [accompanying paper](https://arxiv.org/abs/2412.13043), the ***[tutorial page](https://satfra.github.io/DiFfRG/cpp/TutorialTOC.html)*** in the [documentation](https://satfra.github.io/DiFfRG/cpp/index.html) and the examples in `Examples/`. ## Citation diff --git a/Tutorials/tut1/tut1.cc b/Tutorials/tut1/tut1.cc index 0e13147..ee784ac 100644 --- a/Tutorials/tut1/tut1.cc +++ b/Tutorials/tut1/tut1.cc @@ -30,7 +30,7 @@ int main(int argc, char *argv[]) TimeStepper time_stepper(json, &assembler); // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); // Now we start the timestepping diff --git a/Tutorials/tut2/tut2.cc b/Tutorials/tut2/tut2.cc index 12870ae..24de500 100644 --- a/Tutorials/tut2/tut2.cc +++ b/Tutorials/tut2/tut2.cc @@ -30,7 +30,7 @@ int main(int argc, char *argv[]) TimeStepper time_stepper(json, &assembler); // Set up the initial condition - FlowingVariables initial_condition(discretization); + FE::FlowingVariables initial_condition(discretization); initial_condition.interpolate(model); // Now we start the timestepping diff --git a/VERSION b/VERSION index 3eefcb9..7dea76e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.0.0 +1.0.1 diff --git a/external/build_dealii.sh b/external/build_dealii.sh index 8fc276e..44ce351 100644 --- a/external/build_dealii.sh +++ b/external/build_dealii.sh @@ -14,8 +14,9 @@ cmake -DCMAKE_BUILD_TYPE=DebugRelease \ -DDEAL_II_COMPONENT_DOCUMENTATION=OFF \ -DDEAL_II_ALLOW_PLATFORM_INTROSPECTION=ON \ -DDEAL_II_WITH_CUDA=OFF \ - -DDEAL_II_WITH_MPI=ON \ + -DDEAL_II_WITH_MPI=OFF \ -DDEAL_II_WITH_UMFPACK=ON \ + -DDEAL_II_WITH_TASKFLOW=OFF \ -DDEAL_II_WITH_VTK=OFF \ -DCMAKE_CXX_FLAGS="${CXX_FLAGS}" \ -DCMAKE_EXE_LINKER_FLAGS="${EXE_LINKER_FLAGS}" \ diff --git a/external/build_sundials.sh b/external/build_sundials.sh index d4122e5..be12468 100644 --- a/external/build_sundials.sh +++ b/external/build_sundials.sh @@ -10,7 +10,7 @@ source ./build_scripts/setup_folders.sh cd $BUILD_PATH cmake -DENABLE_OPENMP=ON \ - -DENABLE_MPI=ON \ + -DENABLE_MPI=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_C_FLAGS_RELEASE="${C_FLAGS} -O3 -DNDEBUG" \ -DCMAKE_EXE_LINKER_FLAGS="${EXE_LINKER_FLAGS}" \