diff --git a/DiFfRG/CMakeLists.txt b/DiFfRG/CMakeLists.txt index 52a18b5..e1a2b75 100644 --- a/DiFfRG/CMakeLists.txt +++ b/DiFfRG/CMakeLists.txt @@ -147,4 +147,4 @@ clang_format(format ${headers} ${sources} ${cuda_sources}) file(GLOB_RECURSE CMAKE_FILES CMakeLists.txt) cmake_format(cmake-format ${CMAKE_FILES} - ${CMAKE_CURRENT_SOURCE_DIR}/cmake/setup_build_system.cmake) \ No newline at end of file + ${CMAKE_CURRENT_SOURCE_DIR}/cmake/setup_build_system.cmake) diff --git a/DiFfRG/cmake/setup_build_system.cmake b/DiFfRG/cmake/setup_build_system.cmake index 3b4c317..cc27c1a 100644 --- a/DiFfRG/cmake/setup_build_system.cmake +++ b/DiFfRG/cmake/setup_build_system.cmake @@ -257,16 +257,21 @@ if(USE_CUDA) endif() -cpmaddpackage( - NAME - spdlog - GITHUB_REPOSITORY - gabime/spdlog - VERSION - 1.14.1 - OPTIONS - "CMAKE_BUILD_TYPE Release" - "CMAKE_CXX_FLAGS \"-O3 -DNDEBUG\"") +if(NOT spdlog_ADDED) + cpmaddpackage( + NAME + spdlog + GITHUB_REPOSITORY + gabime/spdlog + VERSION + 1.14.1 + OPTIONS + "CMAKE_BUILD_TYPE Release") + if(spdlog_ADDED) + install(TARGETS spdlog + EXPORT DiFfRGTargets) + endif() +endif() # ############################################################################## # Helper functions diff --git a/DiFfRG/include/DiFfRG/discretization/FV/assembler/KurganovTadmor.hh b/DiFfRG/include/DiFfRG/discretization/FV/assembler/KurganovTadmor.hh new file mode 100644 index 0000000..e0d4c41 --- /dev/null +++ b/DiFfRG/include/DiFfRG/discretization/FV/assembler/KurganovTadmor.hh @@ -0,0 +1,621 @@ +#pragma once + +// external libraries + +// DiFfRG +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +namespace DiFfRG +{ + namespace FV + { + namespace KurganovTadmor + { + using namespace dealii; + + namespace internal + { + /** + * @brief Class to hold data for each assembly thread, i.e. FEValues for cells, interfaces, as well as + * pre-allocated data structures for the solutions + */ + template struct ScratchData { + static constexpr int dim = Discretization::dim; + using NumberType = typename Discretization::NumberType; + using VectorType = Vector; + + ScratchData(const Mapping &mapping, const FiniteElement &fe, const Quadrature &quadrature, + const UpdateFlags update_flags = update_values | update_gradients | update_quadrature_points | + update_JxW_values) + : n_components(fe.n_components()), fe_values(mapping, fe, quadrature, update_flags) + { + solution.resize(quadrature.size(), VectorType(n_components)); + solution_dot.resize(quadrature.size(), VectorType(n_components)); + } + + ScratchData(const ScratchData &scratch_data) + : n_components(scratch_data.fe_values.get_fe().n_components()), + fe_values(scratch_data.fe_values.get_mapping(), scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), scratch_data.fe_values.get_update_flags()) + { + solution.resize(scratch_data.fe_values.get_quadrature().size(), VectorType(n_components)); + solution_dot.resize(scratch_data.fe_values.get_quadrature().size(), VectorType(n_components)); + } + + const uint n_components; + + FEValues fe_values; + + std::vector solution; + std::vector solution_dot; + }; + + // TODO fewer memory allocations + template struct CopyData_R { + struct CopyDataFace_R { + Vector cell_residual; + std::vector joint_dof_indices; + + template void reinit(const FEInterfaceValues &fe_iv) + { + cell_residual.reinit(fe_iv.n_current_interface_dofs()); + joint_dof_indices = fe_iv.get_interface_dof_indices(); + } + }; + + Vector cell_residual; + Vector cell_mass; + std::vector local_dof_indices; + std::vector face_data; + + template void reinit(const Iterator &cell, uint dofs_per_cell) + { + cell_residual.reinit(dofs_per_cell); + cell_mass.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + } + }; + + // TODO fewer memory allocations + template struct CopyData_J { + struct CopyDataFace_J { + FullMatrix cell_jacobian; + FullMatrix extractor_cell_jacobian; + std::vector joint_dof_indices; + + template void reinit(const FEInterfaceValues &fe_iv, uint n_extractors) + { + uint dofs_per_cell = fe_iv.n_current_interface_dofs(); + cell_jacobian.reinit(dofs_per_cell, dofs_per_cell); + if (n_extractors > 0) extractor_cell_jacobian.reinit(dofs_per_cell, n_extractors); + joint_dof_indices = fe_iv.get_interface_dof_indices(); + } + }; + + FullMatrix cell_jacobian; + FullMatrix extractor_cell_jacobian; + FullMatrix cell_mass_jacobian; + std::vector local_dof_indices; + std::vector face_data; + + template void reinit(const Iterator &cell, uint dofs_per_cell, uint n_extractors) + { + cell_jacobian.reinit(dofs_per_cell, dofs_per_cell); + if (n_extractors > 0) extractor_cell_jacobian.reinit(dofs_per_cell, n_extractors); + cell_mass_jacobian.reinit(dofs_per_cell, dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + } + }; + + template struct CopyData_I { + struct CopyFaceData_I { + std::array cell_indices; + std::array values; + }; + std::vector face_data; + double value; + uint cell_index; + }; + } // namespace internal + + template + requires MeshIsRectangular + class Assembler : public AbstractAssembler + { + protected: + template auto fv_tie(T &&...t) + { + return named_tuple, "fe_functions">(std::tie(t...)); + } + + template static constexpr auto v_tie(T &&...t) + { + return named_tuple, "variables", "extractors">(std::tie(t...)); + } + + template static constexpr auto e_tie(T &&...t) + { + return named_tuple, "fe_functions", "fe_derivatives", "fe_hessians", "extractors", + "variables">(std::tie(t...)); + } + + public: + using Discretization = Discretization_; + using Model = Model_; + using NumberType = typename Discretization::NumberType; + using VectorType = typename Discretization::VectorType; + + using Components = typename Discretization::Components; + static constexpr uint dim = Discretization::dim; + static constexpr uint n_components = Components::count_fe_functions(0); + + Assembler(Discretization &discretization, Model &model, const JSONValue &json) + : discretization(discretization), model(model), dof_handler(discretization.get_dof_handler()), + mapping(discretization.get_mapping()), triangulation(discretization.get_triangulation()), json(json), + threads(json.get_uint("/discretization/threads")), + batch_size(json.get_uint("/discretization/batch_size")), + EoM_abs_tol(json.get_double("/discretization/EoM_abs_tol")), + EoM_max_iter(json.get_uint("/discretization/EoM_max_iter")), + quadrature(1 + json.get_uint("/discretization/overintegration")) + { + if (this->threads == 0) this->threads = dealii::MultithreadInfo::n_threads() / 2; + spdlog::get("log")->info("FV: Using {} threads for assembly.", threads); + + reinit(); + } + + virtual void reinit_vector(VectorType &vec) const override + { + const auto block_structure = discretization.get_block_structure(); + vec.reinit(block_structure[0]); + } + + virtual IndexSet get_differential_indices() const override { return IndexSet(); } + + virtual void attach_data_output(DataOutput &data_out, const VectorType &solution, + const VectorType &variables, const VectorType &dt_solution = VectorType(), + const VectorType &residual = VectorType()) + { + const auto fe_function_names = Components::FEFunction_Descriptor::get_names_vector(); + std::vector fe_function_names_residual; + for (const auto &name : fe_function_names) + fe_function_names_residual.push_back(name + "_residual"); + std::vector fe_function_names_dot; + for (const auto &name : fe_function_names) + fe_function_names_dot.push_back(name + "_dot"); + + auto &fe_out = data_out.fe_output(); + fe_out.attach(dof_handler, solution, fe_function_names); + if (dt_solution.size() > 0) fe_out.attach(dof_handler, dt_solution, fe_function_names_dot); + if (residual.size() > 0) fe_out.attach(dof_handler, residual, fe_function_names_residual); + + // readouts(data_out, solution, variables); + } + + virtual void reinit() override + { + Timer timer; + + // Mass sparsity pattern is diagonal + auto dynamic_sparsity_pattern_mass = DynamicSparsityPattern(dof_handler.n_dofs()); + for (uint i = 0; i < dof_handler.n_dofs(); ++i) + dynamic_sparsity_pattern_mass.add(i, i); + sparsity_pattern_mass.copy_from(dynamic_sparsity_pattern_mass); + + mass_matrix = SparseMatrix(sparsity_pattern_mass, IdentityMatrix(dof_handler.n_dofs())); + + constexpr uint stencil = 1; + build_sparsity(sparsity_pattern_jacobian, dof_handler, dof_handler, stencil, true); + + timings_reinit.push_back(timer.wall_time()); + } + + virtual void set_time(double t) override { model.set_time(t); } + + virtual const SparsityPattern &get_sparsity_pattern_jacobian() const override + { + return sparsity_pattern_jacobian; + } + virtual const SparseMatrix &get_mass_matrix() const override { return mass_matrix; } + + virtual void residual_variables(VectorType &residual, const VectorType &variables, const VectorType &) override + { + model.dt_variables(residual, fv_tie(variables)); + }; + + virtual void jacobian_variables(FullMatrix &jacobian, const VectorType &variables, + const VectorType &) override + { + model.template jacobian_variables<0>(jacobian, fv_tie(variables)); + }; + + void readouts(DataOutput &data_out, const VectorType &, const VectorType &variables) const + { + auto helper = [&](auto EoMfun, auto outputter) { + (void)EoMfun; + outputter(data_out, Point<0>(), fv_tie(variables)); + }; + model.template readouts_multiple(helper, data_out); + } + + virtual void mass(VectorType &mass, const VectorType &solution_global, const VectorType &solution_global_dot, + NumberType weight) override + { + using Iterator = typename DoFHandler::active_cell_iterator; + using Scratch = internal::ScratchData; + using CopyData = internal::CopyData_R; + const auto &constraints = discretization.get_constraints(); + + const auto cell_worker = [&](const Iterator &cell, Scratch &scratch_data, CopyData ©_data) { + scratch_data.fe_values.reinit(cell); + const auto &fe_v = scratch_data.fe_values; + const uint n_dofs = fe_v.get_fe().n_dofs_per_cell(); + + copy_data.reinit(cell, n_dofs); + const auto &JxW = fe_v.get_JxW_values(); + const auto &q_points = fe_v.get_quadrature_points(); + const auto &q_indices = fe_v.quadrature_point_indices(); + + auto &solution = scratch_data.solution; + auto &solution_dot = scratch_data.solution_dot; + fe_v.get_function_values(solution_global, solution); + fe_v.get_function_values(solution_global_dot, solution_dot); + + std::array mass{}; + for (const auto &q_index : q_indices) { + const auto &x_q = q_points[q_index]; + model.mass(mass, x_q, solution[q_index], solution_dot[q_index]); + + for (uint i = 0; i < n_dofs; ++i) { + const auto component_i = fe_v.get_fe().system_to_component_index(i).first; + copy_data.cell_residual(i) += weight * JxW[q_index] * + fe_v.shape_value_component(i, q_index, component_i) * + mass[component_i]; // +phi_i(x_q) * mass(x_q, u_q) + } + } + }; + const auto copier = [&](const CopyData &c) { + constraints.distribute_local_to_global(c.cell_residual, c.local_dof_indices, mass); + }; + + Scratch scratch_data(mapping, discretization.get_fe(), quadrature); + CopyData copy_data; + MeshWorker::AssembleFlags flags = MeshWorker::assemble_own_cells; + + MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data, + copy_data, flags, nullptr, nullptr, threads, batch_size); + } + + virtual void residual(VectorType &residual, const VectorType &solution_global, NumberType weight, + const VectorType &solution_global_dot, NumberType weight_mass, + const VectorType &variables = VectorType()) override + { + using Iterator = typename DoFHandler::active_cell_iterator; + using Scratch = internal::ScratchData; + using CopyData = internal::CopyData_R; + const auto &constraints = discretization.get_constraints(); + + const auto cell_worker = [&](const Iterator &cell, Scratch &scratch_data, CopyData ©_data) { + scratch_data.fe_values.reinit(cell); + const auto &fe_v = scratch_data.fe_values; + const uint n_dofs = fe_v.get_fe().n_dofs_per_cell(); + + copy_data.reinit(cell, n_dofs); + const auto &JxW = fe_v.get_JxW_values(); + const auto &q_points = fe_v.get_quadrature_points(); + const auto &q_indices = fe_v.quadrature_point_indices(); + + auto &solution = scratch_data.solution; + auto &solution_dot = scratch_data.solution_dot; + fe_v.get_function_values(solution_global, solution); + fe_v.get_function_values(solution_global_dot, solution_dot); + + std::array mass{}; + std::array source{}; + for (const auto &q_index : q_indices) { + const auto &x_q = q_points[q_index]; + model.mass(mass, x_q, solution[q_index], solution_dot[q_index]); + model.source(source, x_q, fv_tie(solution[q_index])); + + for (uint i = 0; i < n_dofs; ++i) { + const auto component_i = fe_v.get_fe().system_to_component_index(i).first; + copy_data.cell_mass(i) += weight_mass * JxW[q_index] * + fe_v.shape_value_component(i, q_index, component_i) * + mass[component_i]; // +phi_i(x_q) * mass(x_q, u_q) + copy_data.cell_residual(i) += JxW[q_index] * weight * // dx * + (fe_v.shape_value_component(i, q_index, component_i) * + source[component_i]); // -phi_i(x_q) * source(x_q, u_q) + } + } + }; + const auto copier = [&](const CopyData &c) { + constraints.distribute_local_to_global(c.cell_residual, c.local_dof_indices, residual); + constraints.distribute_local_to_global(c.cell_mass, c.local_dof_indices, residual); + }; + + Scratch scratch_data(mapping, discretization.get_fe(), quadrature); + CopyData copy_data; + MeshWorker::AssembleFlags flags = MeshWorker::assemble_own_cells; + + Timer timer; + MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data, + copy_data, flags, nullptr, nullptr, threads, batch_size); + timings_residual.push_back(timer.wall_time()); + } + + virtual void jacobian_mass(SparseMatrix &jacobian, const VectorType &solution_global, + const VectorType &solution_global_dot, NumberType alpha = 1., + NumberType beta = 1.) override + { + using Iterator = typename DoFHandler::active_cell_iterator; + using Scratch = internal::ScratchData; + using CopyData = internal::CopyData_J; + const auto &constraints = discretization.get_constraints(); + + const auto cell_worker = [&](const Iterator &cell, Scratch &scratch_data, CopyData ©_data) { + scratch_data.fe_values.reinit(cell); + const auto &fe_v = scratch_data.fe_values; + const uint n_dofs = fe_v.get_fe().n_dofs_per_cell(); + + copy_data.reinit(cell, n_dofs, Components::count_extractors()); + const auto &JxW = fe_v.get_JxW_values(); + const auto &q_points = fe_v.get_quadrature_points(); + const auto &q_indices = fe_v.quadrature_point_indices(); + + auto &solution = scratch_data.solution; + auto &solution_dot = scratch_data.solution_dot; + fe_v.get_function_values(solution_global, solution); + fe_v.get_function_values(solution_global_dot, solution_dot); + + SimpleMatrix j_mass; + SimpleMatrix j_mass_dot; + for (const auto &q_index : q_indices) { + const auto &x_q = q_points[q_index]; + model.template jacobian_mass<0>(j_mass, x_q, solution[q_index], solution_dot[q_index]); + model.template jacobian_mass<1>(j_mass_dot, x_q, solution[q_index], solution_dot[q_index]); + + for (uint i = 0; i < n_dofs; ++i) { + const auto component_i = fe_v.get_fe().system_to_component_index(i).first; + for (uint j = 0; j < n_dofs; ++j) { + const auto component_j = fe_v.get_fe().system_to_component_index(j).first; + copy_data.cell_jacobian(i, j) += + JxW[q_index] * fe_v.shape_value_component(j, q_index, component_j) * + fe_v.shape_value_component(i, q_index, component_i) * + (alpha * j_mass_dot(component_i, component_j) + beta * j_mass(component_i, component_j)); + } + } + } + }; + const auto copier = [&](const CopyData &c) { + constraints.distribute_local_to_global(c.cell_jacobian, c.local_dof_indices, jacobian); + }; + + Scratch scratch_data(mapping, discretization.get_fe(), quadrature); + CopyData copy_data; + MeshWorker::AssembleFlags flags = MeshWorker::assemble_own_cells; + + Timer timer; + MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data, + copy_data, flags, nullptr, nullptr, threads, batch_size); + timings_jacobian.push_back(timer.wall_time()); + } + + virtual void jacobian(SparseMatrix &jacobian, const VectorType &solution_global, NumberType weight, + const VectorType &solution_global_dot, NumberType alpha, NumberType beta, + const VectorType &variables = VectorType()) override + { + using Iterator = typename DoFHandler::active_cell_iterator; + using Scratch = internal::ScratchData; + using CopyData = internal::CopyData_J; + const auto &constraints = discretization.get_constraints(); + + const auto cell_worker = [&](const Iterator &cell, Scratch &scratch_data, CopyData ©_data) { + scratch_data.fe_values.reinit(cell); + const auto &fe_v = scratch_data.fe_values; + const uint n_dofs = fe_v.get_fe().n_dofs_per_cell(); + + copy_data.reinit(cell, n_dofs, Components::count_extractors()); + const auto &JxW = fe_v.get_JxW_values(); + const auto &q_points = fe_v.get_quadrature_points(); + const auto &q_indices = fe_v.quadrature_point_indices(); + + auto &solution = scratch_data.solution; + auto &solution_dot = scratch_data.solution_dot; + fe_v.get_function_values(solution_global, solution); + fe_v.get_function_values(solution_global_dot, solution_dot); + + SimpleMatrix j_mass; + SimpleMatrix j_mass_dot; + SimpleMatrix j_source; + for (const auto &q_index : q_indices) { + const auto &x_q = q_points[q_index]; + model.template jacobian_mass<0>(j_mass, x_q, solution[q_index], solution_dot[q_index]); + model.template jacobian_mass<1>(j_mass_dot, x_q, solution[q_index], solution_dot[q_index]); + model.template jacobian_source<0, 0>(j_source, x_q, fv_tie(solution[q_index])); + + for (uint i = 0; i < n_dofs; ++i) { + const auto component_i = fe_v.get_fe().system_to_component_index(i).first; + for (uint j = 0; j < n_dofs; ++j) { + const auto component_j = fe_v.get_fe().system_to_component_index(j).first; + copy_data.cell_jacobian(i, j) += + weight * JxW[q_index] * fe_v.shape_value_component(j, q_index, component_j) * // dx * phi_j * ( + (fe_v.shape_value_component(i, q_index, component_i) * + j_source(component_i, component_j)); // -phi_i * jsource) + copy_data.cell_mass_jacobian(i, j) += + JxW[q_index] * fe_v.shape_value_component(j, q_index, component_j) * + fe_v.shape_value_component(i, q_index, component_i) * + (alpha * j_mass_dot(component_i, component_j) + beta * j_mass(component_i, component_j)); + } + } + } + }; + const auto copier = [&](const CopyData &c) { + constraints.distribute_local_to_global(c.cell_jacobian, c.local_dof_indices, jacobian); + constraints.distribute_local_to_global(c.cell_mass_jacobian, c.local_dof_indices, jacobian); + }; + + Scratch scratch_data(mapping, discretization.get_fe(), quadrature); + CopyData copy_data; + MeshWorker::AssembleFlags flags = MeshWorker::assemble_own_cells; + + Timer timer; + MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data, + copy_data, flags, nullptr, nullptr, threads, batch_size); + timings_jacobian.push_back(timer.wall_time()); + } + + void build_sparsity(SparsityPattern &sparsity_pattern, const DoFHandler &to_dofh, + const DoFHandler &from_dofh, const int stencil = 1, + bool add_extractor_dofs = false) const + { + const auto &triangulation = discretization.get_triangulation(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + + const auto to_dofs_per_cell = to_dofh.get_fe().dofs_per_cell; + const auto from_dofs_per_cell = from_dofh.get_fe().dofs_per_cell; + + for (const auto &t_cell : triangulation.active_cell_iterators()) { + std::vector to_dofs(to_dofs_per_cell); + std::vector from_dofs(from_dofs_per_cell); + const auto to_cell = t_cell->as_dof_handler_iterator(to_dofh); + const auto from_cell = t_cell->as_dof_handler_iterator(from_dofh); + to_cell->get_dof_indices(to_dofs); + from_cell->get_dof_indices(from_dofs); + + std::function add_all_neighbor_dofs = [&](const auto &from_cell, + const int stencil_level) { + for (const auto face_no : from_cell->face_indices()) { + const auto face = from_cell->face(face_no); + if (!face->at_boundary()) { + auto neighbor_cell = from_cell->neighbor(face_no); + + if (dim == 1) + while (neighbor_cell->has_children()) + neighbor_cell = neighbor_cell->child(face_no == 0 ? 1 : 0); + + // add all children + else if (neighbor_cell->has_children()) { + throw std::runtime_error("not yet implemented lol"); + } + + if (!neighbor_cell->is_active()) continue; + + std::vector tmp(from_dofs_per_cell); + neighbor_cell->get_dof_indices(tmp); + + from_dofs.insert(std::end(from_dofs), std::begin(tmp), std::end(tmp)); + + if (stencil_level < stencil) add_all_neighbor_dofs(neighbor_cell, stencil_level + 1); + } + } + }; + + add_all_neighbor_dofs(from_cell, 1); + + for (const auto i : to_dofs) + for (const auto j : from_dofs) + dsp.add(i, j); + } + + // if (add_extractor_dofs) + // for (uint row = 0; row < dsp.n_rows(); ++row) + // dsp.add_row_entries(row, extractor_dof_indices); + sparsity_pattern.copy_from(dsp); + } + + void log(const std::string logger) + { + std::stringstream ss; + ss << "FV Assembler: " << std::endl; + ss << " Reinit: " << average_time_reinit() * 1000 << "ms (" << num_reinits() << ")" << std::endl; + ss << " Residual: " << average_time_residual_assembly() * 1000 << "ms (" << num_residuals() << ")" + << std::endl; + ss << " Jacobian: " << average_time_jacobian_assembly() * 1000 << "ms (" << num_jacobians() << ")" + << std::endl; + spdlog::get(logger)->info(ss.str()); + } + + double average_time_reinit() const + { + double t = 0.; + double n = timings_reinit.size(); + for (const auto &t_ : timings_reinit) + t += t_ / n; + return t; + } + uint num_reinits() const { return timings_reinit.size(); } + + double average_time_residual_assembly() const + { + double t = 0.; + double n = timings_residual.size(); + for (const auto &t_ : timings_residual) + t += t_ / n; + return t; + } + uint num_residuals() const { return timings_residual.size(); } + + double average_time_jacobian_assembly() const + { + double t = 0.; + double n = timings_jacobian.size(); + for (const auto &t_ : timings_jacobian) + t += t_ / n; + return t; + } + uint num_jacobians() const { return timings_jacobian.size(); } + + protected: + Discretization &discretization; + Model &model; + const DoFHandler &dof_handler; + const Mapping &mapping; + const Triangulation &triangulation; + + const JSONValue &json; + + uint threads; + const uint batch_size; + + mutable typename Triangulation::cell_iterator EoM_cell; + typename Triangulation::cell_iterator old_EoM_cell; + const double EoM_abs_tol; + const uint EoM_max_iter; + + const QGauss quadrature; + + SparsityPattern sparsity_pattern_mass; + SparsityPattern sparsity_pattern_jacobian; + SparseMatrix mass_matrix; + + std::vector timings_reinit; + std::vector timings_residual; + std::vector timings_jacobian; + }; + } // namespace KurganovTadmor + } // namespace FV +} // namespace DiFfRG \ No newline at end of file diff --git a/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh b/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh index 741f39a..77a3967 100644 --- a/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh +++ b/DiFfRG/include/DiFfRG/discretization/FV/discretization.hh @@ -5,14 +5,14 @@ #include #include #include -#include +#include #include #include #include -#include // DiFfRG #include +#include namespace DiFfRG { @@ -22,7 +22,7 @@ namespace DiFfRG /** * @brief Class to manage the system on which we solve, i.e. fe spaces, grids, etc. - * This class is a System for CG systems. + * This class is a System for FV systems. * * @tparam Model_ The Model class used for the Simulation */ @@ -36,21 +36,93 @@ namespace DiFfRG using Mesh = Mesh_; static constexpr uint dim = Mesh::dim; - Discretization(Mesh &mesh, const JSONValue &json) : mesh(mesh), json(json) {}; + Discretization(Mesh &mesh, const JSONValue &json) + : mesh(mesh), json(json), + fe(std::make_shared>(FE_DGQ(0), Components::count_fe_functions(0))), + dof_handler(mesh.get_triangulation()) + { + setup_dofs(); + } + const auto &get_constraints(const uint i = 0) const + { + (void)i; + return constraints; + } + auto &get_constraints(const uint i = 0) + { + (void)i; + return constraints; + } + const auto &get_dof_handler(const uint i = 0) const + { + (void)i; + return dof_handler; + } + auto &get_dof_handler(const uint i = 0) + { + (void)i; + return dof_handler; + } + const auto &get_fe(uint i = 0) const + { + if (i != 0) throw std::runtime_error("Wrong FE index"); + return *fe; + } const auto &get_mapping() const { return mapping; } const auto &get_triangulation() const { return mesh.get_triangulation(); } auto &get_triangulation() { return mesh.get_triangulation(); } + const Point &get_support_point(const uint &dof) const { return support_points[dof]; } + const auto &get_support_points() const { return support_points; } const auto &get_json() const { return json; } - void reinit() {} + void reinit() { setup_dofs(); } + + uint get_closest_dof(const Point &p) const + { + uint dof = 0; + double min_dist = std::numeric_limits::max(); + for (uint i = 0; i < support_points.size(); ++i) { + const auto dist = p.distance(support_points[i]); + if (dist < min_dist) { + min_dist = dist; + dof = i; + } + } + 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() + { + dof_handler.distribute_dofs(*fe); + + spdlog::get("log")->info("FV: Number of active cells: {}", mesh.get_triangulation().n_active_cells()); + spdlog::get("log")->info("FV: Number of degrees of freedom: {}", dof_handler.n_dofs()); + + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + support_points.resize(dof_handler.n_dofs()); + DoFTools::map_dofs_to_support_points(mapping, dof_handler, support_points); + } + Mesh &mesh; JSONValue json; + std::shared_ptr> fe; + DoFHandler dof_handler; AffineConstraints constraints; MappingQ1 mapping; + std::vector> support_points; }; } // namespace FV -} // namespace DiFfRG +} // namespace DiFfRG \ No newline at end of file diff --git a/DiFfRG/include/DiFfRG/discretization/common/types.hh b/DiFfRG/include/DiFfRG/discretization/common/types.hh new file mode 100644 index 0000000..effd2c5 --- /dev/null +++ b/DiFfRG/include/DiFfRG/discretization/common/types.hh @@ -0,0 +1,7 @@ +#pragma once + +namespace DiFfRG +{ + template + concept MeshIsRectangular = requires(T x) { T::is_rectangular; } && T::is_rectangular; +} // namespace DiFfRG \ No newline at end of file diff --git a/DiFfRG/include/DiFfRG/discretization/data/data.hh b/DiFfRG/include/DiFfRG/discretization/data/data.hh index 108dee2..50264c0 100644 --- a/DiFfRG/include/DiFfRG/discretization/data/data.hh +++ b/DiFfRG/include/DiFfRG/discretization/data/data.hh @@ -168,73 +168,8 @@ namespace DiFfRG 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; - }; + using FlowingVariables = DiFfRG::FE::FlowingVariables; } // namespace FV } // namespace DiFfRG diff --git a/DiFfRG/include/DiFfRG/discretization/discretization.hh b/DiFfRG/include/DiFfRG/discretization/discretization.hh index d20bdeb..04bd5d6 100644 --- a/DiFfRG/include/DiFfRG/discretization/discretization.hh +++ b/DiFfRG/include/DiFfRG/discretization/discretization.hh @@ -1,5 +1,7 @@ #pragma once +#include + #include #include #include diff --git a/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh b/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh index 0e74676..c8681a9 100644 --- a/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh +++ b/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh @@ -20,6 +20,7 @@ namespace DiFfRG { public: static constexpr uint dim = dim_; + static constexpr bool is_rectangular = true; /** * @brief Construct a new RectangularMesh object. diff --git a/DiFfRG/src/CMakeLists.txt b/DiFfRG/src/CMakeLists.txt index 475531b..727164d 100644 --- a/DiFfRG/src/CMakeLists.txt +++ b/DiFfRG/src/CMakeLists.txt @@ -1,6 +1,5 @@ # add_compile_definitions(CUDA_API_PER_THREAD_DEFAULT_STREAM) -add_library( - DiFfRG STATIC +set(DiFfRG_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/common/boost_json.cc ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/common/configuration_helper.cc ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/common/csv_reader.cc @@ -30,9 +29,7 @@ add_library( ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/timestepping/trbdf2.cc ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/flow_equations.cc) if(USE_CUDA) - target_sources( - DiFfRG - PRIVATE + list(APPEND DiFfRG_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/common/cuda_prefix.cu ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/integration/quadrature_provider.cu ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/interpolation/linear_interpolation_1D.cu @@ -43,9 +40,7 @@ if(USE_CUDA) ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/interpolation/tex_linear_interpolation_1D_stack.cu ) else() - target_sources( - DiFfRG - PRIVATE + list(APPEND DiFfRG_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/integration/quadrature_provider.cc ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/interpolation/linear_interpolation_1D.cc ${CMAKE_CURRENT_SOURCE_DIR}/DiFfRG/physics/interpolation/linear_interpolation_2D.cc @@ -56,6 +51,8 @@ else() ) endif() +add_library(DiFfRG STATIC ${DiFfRG_SOURCES}) + setup_target(DiFfRG) target_include_directories( diff --git a/build.sh b/build.sh index d2c3753..ab5778d 100644 --- a/build.sh +++ b/build.sh @@ -108,7 +108,7 @@ if [[ ${install_dir} != "n" ]] && [[ ${install_dir} != "N" ]]; then # Make sure the install directory is absolute idir=$(expandPath ${install_dir}/) - idir=$(readlink --canonicalize ${idir}) + #idir=$(readlink --canonicalize ${idir}) echo "DiFfRG library will be installed in ${idir}" # Check if the install directory is writable diff --git a/external/build_boost.sh b/external/build_boost.sh index 6e77b5b..82c3856 100644 --- a/external/build_boost.sh +++ b/external/build_boost.sh @@ -28,8 +28,8 @@ fi echo $COMPILER cd $SOURCE_PATH -./bootstrap.sh --prefix=${INSTALL_PATH} -./b2 --build-dir=${BUILD_PATH} \ +$SuperUser ./bootstrap.sh --prefix=${INSTALL_PATH} +$SuperUser ./b2 --build-dir=${BUILD_PATH} \ --prefix=${INSTALL_PATH} \ --with-headers \ --with-iostreams \ diff --git a/external/build_dealii.sh b/external/build_dealii.sh index 45fbf74..9b11da7 100644 --- a/external/build_dealii.sh +++ b/external/build_dealii.sh @@ -33,4 +33,4 @@ cmake -DCMAKE_BUILD_TYPE=DebugRelease \ 2>&1 | tee $CMAKE_LOG_FILE make -j $THREADS 2>&1 | tee $MAKE_LOG_FILE -make -j $THREADS install +$SuperUser make -j $THREADS install diff --git a/external/build_kokkos.sh b/external/build_kokkos.sh index bded9fd..e810945 100644 --- a/external/build_kokkos.sh +++ b/external/build_kokkos.sh @@ -13,15 +13,15 @@ source $SCRIPT_PATH/build_scripts/setup_folders.sh cd $BUILD_PATH - #-DKokkos_ARCH_NATIVE=ON \ +#-DKokkos_ARCH_NATIVE=ON \ cmake -DKokkos_ENABLE_SERIAL=ON \ - -DKokkos_ENABLE_OPENMP=OFF \ - -DKokkos_ENABLE_CUDA=OFF \ - -DCMAKE_CXX_FLAGS="${CXX_FLAGS}" \ - -DCMAKE_EXE_LINKER_FLAGS="${EXE_LINKER_FLAGS}" \ - -DCMAKE_INSTALL_PREFIX=${INSTALL_PATH} \ - -S ${SOURCE_PATH} \ - 2>&1 | tee $CMAKE_LOG_FILE + -DKokkos_ENABLE_OPENMP=OFF \ + -DKokkos_ENABLE_CUDA=OFF \ + -DCMAKE_CXX_FLAGS="${CXX_FLAGS}" \ + -DCMAKE_EXE_LINKER_FLAGS="${EXE_LINKER_FLAGS}" \ + -DCMAKE_INSTALL_PREFIX=${INSTALL_PATH} \ + -S ${SOURCE_PATH} \ + 2>&1 | tee $CMAKE_LOG_FILE make -j $THREADS 2>&1 | tee $MAKE_LOG_FILE -make -j $THREADS install +$SuperUser make -j $THREADS install diff --git a/external/build_oneTBB.sh b/external/build_oneTBB.sh index c1771a3..38d971a 100644 --- a/external/build_oneTBB.sh +++ b/external/build_oneTBB.sh @@ -23,4 +23,4 @@ cmake -DCMAKE_BUILD_TYPE=Release \ 2>&1 | tee $CMAKE_LOG_FILE make -j $THREADS 2>&1 | tee $MAKE_LOG_FILE -make -j $THREADS install +$SuperUser make -j $THREADS install diff --git a/external/build_scripts/setup_folders.sh b/external/build_scripts/setup_folders.sh index 9b2c685..1c2e7c6 100644 --- a/external/build_scripts/setup_folders.sh +++ b/external/build_scripts/setup_folders.sh @@ -1,3 +1,17 @@ mkdir -p $SCRIPT_PATH/logs mkdir -p $BUILD_PATH -mkdir -p $INSTALL_PATH + +failed_first=0 +mkdir -p $INSTALL_PATH &>/dev/null && touch $INSTALL_PATH/_permission_test &>/dev/null || { + failed_first=1 +} + +SuperUser="" +if ((failed_first == 0)); then + echo "Elevated permissions required to write into $INSTALL_PATH." + SuperUser="sudo -E" + + $SuperUser mkdir -p $INSTALL_PATH +else + rm $INSTALL_PATH/_permission_test +fi diff --git a/external/build_sundials.sh b/external/build_sundials.sh index 1ac76eb..d77ef0c 100644 --- a/external/build_sundials.sh +++ b/external/build_sundials.sh @@ -23,4 +23,4 @@ cmake -DENABLE_OPENMP=ON \ 2>&1 | tee $CMAKE_LOG_FILE make -j $THREADS 2>&1 | tee $MAKE_LOG_FILE -make -j $THREADS install +$SuperUser make -j $THREADS install diff --git a/update_DiFfRG.sh b/update_DiFfRG.sh index c632e67..3d1e96c 100644 --- a/update_DiFfRG.sh +++ b/update_DiFfRG.sh @@ -123,7 +123,7 @@ if [[ ${option_install_library} != "n" ]] && [[ ${option_install_library} != "N" # Make sure the install directory is absolute idir=$(expandPath ${option_install_library}/) - idir=$(readlink --canonicalize ${idir}) + #idir=$(readlink --canonicalize ${idir}) echo "DiFfRG library will be installed in ${idir}" # Check if the install directory is writable