Skip to content
Snippets Groups Projects
Commit 3a2cafbb authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Suppress Connectivity1D

- Connectivity1D is now Connectivity<1>

- Connectivity1D can no longer be built using a number of cells (this will be
   possible again with the DSL)
parent d6369864
No related branches found
No related tags found
1 merge request!6Feature/crs
...@@ -11,7 +11,6 @@ ...@@ -11,7 +11,6 @@
// #include <AcousticSolverClass.hpp> // #include <AcousticSolverClass.hpp>
// #include <AcousticSolverTest.hpp> // #include <AcousticSolverTest.hpp>
#include <Connectivity1D.hpp>
#include <Connectivity3D.hpp> #include <Connectivity3D.hpp>
#include <Mesh.hpp> #include <Mesh.hpp>
...@@ -441,89 +440,8 @@ int main(int argc, char *argv[]) ...@@ -441,89 +440,8 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n"; std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
} else { } else {
// class for acoustic solver test std::cerr << "Connectivity1D defined by number of nodes no more implemented\n";
Kokkos::Timer timer; std::exit(0);
timer.reset();
std::shared_ptr<Connectivity1D>connectivity( new Connectivity1D(number));
typedef Mesh<Connectivity1D> MeshType;
typedef MeshData<MeshType> MeshDataType;
typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
MeshType mesh(connectivity);
MeshDataType mesh_data(mesh);
std::vector<BoundaryConditionHandler> bc_list;
{ // quite dirty!
// SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0
// = new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}),
// TinyVector<1>(-1));
// std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0);
// bc_list.push_back(BoundaryConditionHandler(bc0));
PressureBoundaryCondition* pres_bc0
= new PressureBoundaryCondition(1,
std::vector<unsigned int>({static_cast<unsigned int>(0u)}));
std::shared_ptr<PressureBoundaryCondition> bc0(pres_bc0);
bc_list.push_back(BoundaryConditionHandler(bc0));
PressureBoundaryCondition* pres_bc1
= new PressureBoundaryCondition(0.1,
std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
bc_list.push_back(BoundaryConditionHandler(bc1));
}
UnknownsType unknowns(mesh_data);
unknowns.initializeSod();
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
typedef TinyVector<MeshType::dimension> Rd;
const Kokkos::View<const double*> Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
int itermax=std::numeric_limits<int>::max();
int iteration=0;
Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<double*> ej = unknowns.ej();
Kokkos::View<double*> pj = unknowns.pj();
Kokkos::View<double*> gammaj = unknowns.gammaj();
Kokkos::View<double*> cj = unknowns.cj();
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
while((t<tmax) and (iteration<itermax)) {
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
}
acoustic_solver.computeNextStep(t,dt, unknowns);
block_eos.updatePandCFromRhoE();
t += dt;
++iteration;
}
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
{ // gnuplot output for density
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj();
std::ofstream fout("rho");
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << rhoj[j] << '\n';
}
}
} }
} }
catch (const AssertError& error) { catch (const AssertError& error) {
......
#ifndef CONNECTIVITY_1D_HPP
#define CONNECTIVITY_1D_HPP
#include <Kokkos_Core.hpp>
#include <PastisAssert.hpp>
#include <TinyVector.hpp>
#include <ConnectivityUtils.hpp>
#include <TypeOfItem.hpp>
#include <RefId.hpp>
#include <RefNodeList.hpp>
class Connectivity1D
{
public:
static constexpr size_t dimension = 1;
ConnectivityMatrix m_cell_to_node_matrix;
ConnectivityMatrix m_node_to_cell_matrix;
ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
// Stores numbering of nodes of each cell.
//
// This is different from m_cell_to_node_matrix which return the global id of
// a local node in a cell
ConnectivityMatrix m_node_id_per_cell_matrix;
ConnectivityMatrix subItemIdPerItemMatrix() const
{
return m_node_id_per_cell_matrix;
}
private:
std::vector<RefNodeList> m_ref_node_list;
Kokkos::View<double*> m_inv_cell_nb_nodes;
std::vector<std::vector<unsigned int>>
_buildConnectivity(const size_t& number_of_cells)
{
std::vector<std::vector<unsigned int>> cell_by_node_vector(number_of_cells);
for (unsigned int j=0; j<number_of_cells; ++j) {
cell_by_node_vector[j] = {j, j+1};
}
return cell_by_node_vector;
}
public:
void addRefNodeList(const RefNodeList& ref_node_list)
{
m_ref_node_list.push_back(ref_node_list);
}
size_t numberOfRefNodeList() const
{
return m_ref_node_list.size();
}
const RefNodeList& refNodeList(const size_t& i) const
{
return m_ref_node_list[i];
}
KOKKOS_INLINE_FUNCTION
size_t numberOfNodes() const
{
return m_node_to_cell_matrix.numRows();
}
KOKKOS_INLINE_FUNCTION
size_t numberOfCells() const
{
return m_cell_to_node_matrix.numRows();
}
const Kokkos::View<const double*> invCellNbNodes() const
{
return m_inv_cell_nb_nodes;
}
Connectivity1D(const Connectivity1D&) = delete;
Connectivity1D(const size_t& number_of_cells)
: Connectivity1D(_buildConnectivity(number_of_cells))
{
;
}
Connectivity1D(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
{
m_cell_to_node_matrix
= Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
cell_by_node_vector);
Assert(this->numberOfCells()>0);
{
Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const size_t& j) {
const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
inv_cell_nb_nodes[j] = 1./cell_nodes.length;
});
m_inv_cell_nb_nodes = inv_cell_nb_nodes;
}
{
std::vector<std::vector<unsigned int>> node_id_per_cell_vector(this->numberOfCells());
unsigned int id=0;
for (unsigned int j=0; j<this->numberOfCells(); ++j) {
const auto& cell_to_node = m_cell_to_node_matrix.rowConst(j);
auto& node_id_per_cell = node_id_per_cell_vector[j];
node_id_per_cell.resize(cell_to_node.length);
for (size_t r=0; r<cell_to_node.length; ++r) {
node_id_per_cell[r] = id++;
}
}
m_node_id_per_cell_matrix
= Kokkos::create_staticcrsgraph<ConnectivityMatrix>("node_id_per_cell_matrix",
cell_by_node_vector);
}
ConnectivityUtils utils;
utils.computeNodeCellConnectivity(m_cell_to_node_matrix,
m_node_to_cell_matrix,
m_node_to_cell_local_node_matrix);
}
~Connectivity1D()
{
;
}
};
#endif // CONNECTIVITY_1D_HPP
...@@ -24,6 +24,18 @@ class Connectivity; ...@@ -24,6 +24,18 @@ class Connectivity;
template <size_t Dimension> template <size_t Dimension>
class ConnectivityFace; class ConnectivityFace;
template<>
class ConnectivityFace<1>
{
public:
friend struct Hash;
struct Hash
{
size_t operator()(const ConnectivityFace& f) const;
};
};
template<> template<>
class ConnectivityFace<2> class ConnectivityFace<2>
{ {
...@@ -364,8 +376,10 @@ private: ...@@ -364,8 +376,10 @@ private:
m_node_to_cell_matrix, m_node_to_cell_matrix,
m_node_to_cell_local_node_matrix); m_node_to_cell_local_node_matrix);
if constexpr (Dimension>1) {
this->_computeFaceCellConnectivities(); this->_computeFaceCellConnectivities();
} }
}
~Connectivity() ~Connectivity()
{ {
...@@ -663,4 +677,15 @@ Connectivity<2>::subItemIdPerItemMatrix<TypeOfItem::node, ...@@ -663,4 +677,15 @@ Connectivity<2>::subItemIdPerItemMatrix<TypeOfItem::node,
return m_node_id_per_cell_matrix; return m_node_id_per_cell_matrix;
} }
using Connectivity1D = Connectivity<1>;
template <>
template <>
inline const ConnectivityMatrix&
Connectivity<1>::subItemIdPerItemMatrix<TypeOfItem::node,
TypeOfItem::cell>() const
{
return m_node_id_per_cell_matrix;
}
#endif // CONNECTIVITY_3D_HPP #endif // CONNECTIVITY_3D_HPP
...@@ -5,7 +5,6 @@ ...@@ -5,7 +5,6 @@
#include <set> #include <set>
#include <rang.hpp> #include <rang.hpp>
#include <Connectivity1D.hpp>
#include <Connectivity3D.hpp> #include <Connectivity3D.hpp>
#include <Mesh.hpp> #include <Mesh.hpp>
......
...@@ -164,11 +164,7 @@ class SubItemValuePerItem ...@@ -164,11 +164,7 @@ class SubItemValuePerItem
template <typename ConnectivityType> template <typename ConnectivityType>
SubItemValuePerItem(const ConnectivityType& connectivity) SubItemValuePerItem(const ConnectivityType& connectivity)
{ {
if constexpr (ConnectivityType::dimension > 1) {
m_subitem_id_per_item_matrix = connectivity.template subItemIdPerItemMatrix<SubItemType,ItemType>(); m_subitem_id_per_item_matrix = connectivity.template subItemIdPerItemMatrix<SubItemType,ItemType>();
} else {
m_subitem_id_per_item_matrix = connectivity.subItemIdPerItemMatrix();
}
m_values = Kokkos::View<DataType*>("values", m_subitem_id_per_item_matrix.entries.extent(0)); m_values = Kokkos::View<DataType*>("values", m_subitem_id_per_item_matrix.entries.extent(0));
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment