diff --git a/src/main.cpp b/src/main.cpp
index fbbdfd31e136c4318c9463382a245a6b714e1770..4c7c499650008c33f8118ad42e1ab84e98653224 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -8,6 +8,7 @@
#include <Mesh.hpp>
#include <BoundaryCondition.hpp>
#include <AcousticSolver.hpp>
+#include <VoronoiSolver.hpp>
#include <VTKWriter.hpp>
@@ -41,8 +42,19 @@ int main(int argc, char *argv[])
std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
- switch (p_mesh->meshDimension()) {
- case 1: {
+ const bool use_voronoi
+ = [&] () {
+ char* voronoi_env = getenv("VORONOI");
+ return voronoi_env != nullptr;
+ } ();
+
+ if (use_voronoi) {
+ if (p_mesh->meshDimension() != 1) {
+ std::cerr << "higher dimension than 1 is not implemented for Voronoi\n";
+ std::abort();
+ }
+
+ {
std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
for (const auto& sym_boundary_name : sym_boundary_name_list){
@@ -97,7 +109,7 @@ int main(int argc, char *argv[])
unknowns.initializeSod();
- AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
+ VoronoiSolver<MeshDataType> voronoi_solver(mesh_data, unknowns, bc_list);
using Rd = TinyVector<MeshType::dimension>;
@@ -120,13 +132,14 @@ int main(int argc, char *argv[])
VTKWriter vtk_writer("mesh", 0.01);
while((t<tmax) and (iteration<itermax)) {
+ std::cout << "time " << t << ' ';
vtk_writer.write(mesh, t);
- double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt = 0.4*voronoi_solver.voronoi_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
}
- acoustic_solver.computeNextStep(t,dt, unknowns);
-
+ voronoi_solver.computeNextStep(t,dt, unknowns);
+ std::cout << "dt=" << dt << '\n';
block_eos.updatePandCFromRhoE();
t += dt;
@@ -135,9 +148,9 @@ int main(int argc, char *argv[])
vtk_writer.write(mesh, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+ << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+ method_cost_map["VoronoiSolverWithMesh"] = timer.seconds();
{ // gnuplot output for density
const CellValue<const Rd>& xj = mesh_data.xj();
@@ -148,202 +161,312 @@ int main(int argc, char *argv[])
}
}
- break;
- }
- case 2: {
- // test case boundary condition description
- std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
- std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
- for (const auto& sym_boundary_name : sym_boundary_name_list){
- std::shared_ptr<BoundaryDescriptor> boudary_descriptor
- = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
- SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
- = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
- bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
- }
-
- using ConnectivityType = Connectivity2D;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<MeshType>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
-
- const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
- Timer timer;
- timer.reset();
- MeshDataType mesh_data(mesh);
+ }
+ } else {
+ switch (p_mesh->meshDimension()) {
+ case 1: {
+ std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
+ std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+ for (const auto& sym_boundary_name : sym_boundary_name_list){
+ std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+ = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+ SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+ = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+ bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+ }
- std::vector<BoundaryConditionHandler> bc_list;
- {
- for (const auto& bc_descriptor : bc_descriptor_list) {
- switch (bc_descriptor->type()) {
- case BoundaryConditionDescriptor::Type::symmetry: {
- const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
- = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
- for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
- ++i_ref_face_list) {
- const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == sym_bc_descriptor.boundaryDescriptor()) {
- SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
- = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list));
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
+ using ConnectivityType = Connectivity1D;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<MeshType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+ Timer timer;
+ timer.reset();
+ MeshDataType mesh_data(mesh);
+
+ std::vector<BoundaryConditionHandler> bc_list;
+ {
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ switch (bc_descriptor->type()) {
+ case BoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
+ = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList();
+ ++i_ref_node_list) {
+ const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list);
+ const RefId& ref = ref_node_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+ = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_node_list));
+ std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+ bc_list.push_back(BoundaryConditionHandler(bc));
+ }
}
+ break;
+ }
+ default: {
+ perr() << "Unknown BCDescription\n";
+ std::exit(1);
}
- break;
- }
- default: {
- perr() << "Unknown BCDescription\n";
- std::exit(1);
}
}
}
- }
- UnknownsType unknowns(mesh_data);
+ UnknownsType unknowns(mesh_data);
- unknowns.initializeSod();
+ unknowns.initializeSod();
- AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
+ AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const CellValue<const double>& Vj = mesh_data.Vj();
+ using Rd = TinyVector<MeshType::dimension>;
- const double tmax=0.2;
- double t=0;
+ const CellValue<const double>& Vj = mesh_data.Vj();
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
+ const double tmax=0.2;
+ double t=0;
- CellValue<double>& rhoj = unknowns.rhoj();
- CellValue<double>& ej = unknowns.ej();
- CellValue<double>& pj = unknowns.pj();
- CellValue<double>& gammaj = unknowns.gammaj();
- CellValue<double>& cj = unknowns.cj();
+ int itermax=std::numeric_limits<int>::max();
+ int iteration=0;
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
- VTKWriter vtk_writer("mesh", 0.01);
+ BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
- while((t<tmax) and (iteration<itermax)) {
- vtk_writer.write(mesh, t);
- double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
- if (t+dt>tmax) {
- dt=tmax-t;
+ VTKWriter vtk_writer("mesh", 0.01);
+
+ while((t<tmax) and (iteration<itermax)) {
+ vtk_writer.write(mesh, t);
+ 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;
}
- acoustic_solver.computeNextStep(t,dt, unknowns);
+ vtk_writer.write(mesh, t, true); // forces last output
- block_eos.updatePandCFromRhoE();
+ pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
+ << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
- t += dt;
- ++iteration;
+ method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+
+ { // gnuplot output for density
+ const CellValue<const Rd>& xj = mesh_data.xj();
+ const CellValue<const double>& rhoj = unknowns.rhoj();
+ std::ofstream fout("rho");
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+ }
+ }
+
+ break;
}
- vtk_writer.write(mesh, t, true); // forces last output
+ case 2: {
+ // test case boundary condition description
+ std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
+ std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+ for (const auto& sym_boundary_name : sym_boundary_name_list){
+ std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+ = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+ SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+ = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+ bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+ }
- pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+ using ConnectivityType = Connectivity2D;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<MeshType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+ Timer timer;
+ timer.reset();
+ MeshDataType mesh_data(mesh);
+
+ std::vector<BoundaryConditionHandler> bc_list;
+ {
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ switch (bc_descriptor->type()) {
+ case BoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
+ = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+ ++i_ref_face_list) {
+ const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+ = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list));
+ std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+ bc_list.push_back(BoundaryConditionHandler(bc));
+ }
+ }
+ break;
+ }
+ default: {
+ perr() << "Unknown BCDescription\n";
+ std::exit(1);
+ }
+ }
+ }
+ }
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
- break;
- }
- case 3: {
- std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
- std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
- for (const auto& sym_boundary_name : sym_boundary_name_list){
- std::shared_ptr<BoundaryDescriptor> boudary_descriptor
- = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
- SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
- = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+ UnknownsType unknowns(mesh_data);
- bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
- }
+ unknowns.initializeSod();
- using ConnectivityType = Connectivity3D;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<MeshType>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+ AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+ const CellValue<const double>& Vj = mesh_data.Vj();
- Timer timer;
- timer.reset();
- MeshDataType mesh_data(mesh);
+ const double tmax=0.2;
+ double t=0;
- std::vector<BoundaryConditionHandler> bc_list;
- {
- for (const auto& bc_descriptor : bc_descriptor_list) {
- switch (bc_descriptor->type()) {
- case BoundaryConditionDescriptor::Type::symmetry: {
- const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
- = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
- for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
- ++i_ref_face_list) {
- const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == sym_bc_descriptor.boundaryDescriptor()) {
- SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
- = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list));
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
+ int itermax=std::numeric_limits<int>::max();
+ int iteration=0;
+
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
+
+ BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+
+ VTKWriter vtk_writer("mesh", 0.01);
+
+ while((t<tmax) and (iteration<itermax)) {
+ vtk_writer.write(mesh, t);
+ 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;
+ }
+ vtk_writer.write(mesh, t, true); // forces last output
+
+ pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
+ << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+ method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+ break;
+ }
+ case 3: {
+ std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+ std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+ for (const auto& sym_boundary_name : sym_boundary_name_list){
+ std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+ = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+ SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+ = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+ bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+ }
+
+ using ConnectivityType = Connectivity3D;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<MeshType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+ Timer timer;
+ timer.reset();
+ MeshDataType mesh_data(mesh);
+
+ std::vector<BoundaryConditionHandler> bc_list;
+ {
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ switch (bc_descriptor->type()) {
+ case BoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
+ = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+ ++i_ref_face_list) {
+ const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+ = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list));
+ std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+ bc_list.push_back(BoundaryConditionHandler(bc));
+ }
}
+ break;
+ }
+ default: {
+ perr() << "Unknown BCDescription\n";
+ std::exit(1);
}
- break;
- }
- default: {
- perr() << "Unknown BCDescription\n";
- std::exit(1);
}
}
}
- }
- UnknownsType unknowns(mesh_data);
+ UnknownsType unknowns(mesh_data);
- unknowns.initializeSod();
+ unknowns.initializeSod();
- AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
+ AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const CellValue<const double>& Vj = mesh_data.Vj();
+ const CellValue<const double>& Vj = mesh_data.Vj();
- const double tmax=0.2;
- double t=0;
+ const double tmax=0.2;
+ double t=0;
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
+ int itermax=std::numeric_limits<int>::max();
+ int iteration=0;
- CellValue<double>& rhoj = unknowns.rhoj();
- CellValue<double>& ej = unknowns.ej();
- CellValue<double>& pj = unknowns.pj();
- CellValue<double>& gammaj = unknowns.gammaj();
- CellValue<double>& cj = unknowns.cj();
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+ BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
- VTKWriter vtk_writer("mesh", 0.01);
+ VTKWriter vtk_writer("mesh", 0.01);
- while((t<tmax) and (iteration<itermax)) {
- vtk_writer.write(mesh, t);
- 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();
+ while((t<tmax) and (iteration<itermax)) {
+ vtk_writer.write(mesh, t);
+ 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;
- }
- vtk_writer.write(mesh, t, true); // forces last output
+ t += dt;
+ ++iteration;
+ }
+ vtk_writer.write(mesh, t, true); // forces last output
- pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+ pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
+ << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
- break;
+ method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+ break;
+ }
}
}
-
pout() << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
} else {
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 6bd7e892576d329a82de8e10519f5b29fabb6e57..f9a2aeb575ac5265073d47ce2c4025da39db7d51 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -179,7 +179,7 @@ class MeshData
const auto& face_nodes = face_to_node_matrix[l];
#warning should this lambda be replaced by a precomputed correspondance?
- std::function local_node_number_in_cell
+ auto local_node_number_in_cell
= [&](const NodeId& node_number) {
for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
if (node_number == cell_nodes[i_node]) {
diff --git a/src/scheme/VoronoiSolver.hpp b/src/scheme/VoronoiSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..831ed235b70535e0ba279b4413ad2cea4d19f7ea
--- /dev/null
+++ b/src/scheme/VoronoiSolver.hpp
@@ -0,0 +1,392 @@
+#ifndef VORONOI_SOLVER_HPP
+#define VORONOI_SOLVER_HPP
+
+#include <rang.hpp>
+
+#include <ArrayUtils.hpp>
+
+#include <BlockPerfectGas.hpp>
+#include <PastisAssert.hpp>
+
+#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+#include <Mesh.hpp>
+#include <MeshData.hpp>
+#include <FiniteVolumesEulerUnknowns.hpp>
+#include <BoundaryCondition.hpp>
+
+#include <SubItemValuePerItem.hpp>
+
+template<typename MeshData>
+class VoronoiSolver
+{
+ using MeshType = typename MeshData::MeshType;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshData>;
+
+ MeshData& m_mesh_data;
+ const MeshType& m_mesh;
+ const typename MeshType::Connectivity& m_connectivity;
+ const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
+
+ constexpr static size_t dimension = MeshType::dimension;
+
+ using Rd = TinyVector<dimension>;
+ using Rdd = TinyMatrix<dimension>;
+
+ private:
+ PASTIS_INLINE
+ const CellValue<const double>
+ computeRhoCj(const CellValue<const double>& rhoj,
+ const CellValue<const double>& cj)
+ {
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ m_rhocj[j] = rhoj[j]*cj[j];
+ });
+ return m_rhocj;
+ }
+
+ PASTIS_INLINE
+ void computeAjr(const CellValue<const double>& rhocj,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const NodeValuePerCell<const double>& ljr,
+ const NodeValuePerCell<const Rd>& njr)
+ {
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
+ const double& rho_c = rhocj[j];
+ for (size_t r=0; r<nb_nodes; ++r) {
+ m_Ajr(j,r) = tensorProduct(rho_c*Cjr(j,r), njr(j,r));
+ }
+ });
+ }
+
+ PASTIS_INLINE
+ const NodeValue<const Rdd>
+ computeAr(const NodeValuePerCell<const Rdd>& Ajr) {
+ const auto& node_to_cell_matrix
+ = m_connectivity.nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells
+ = m_connectivity.nodeLocalNumbersInTheirCells();
+
+ parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+ Rdd sum = zero;
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells
+ = node_local_numbers_in_their_cells.itemValues(r);
+
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ sum += Ajr(J,R);
+ }
+ m_Ar[r] = sum;
+ });
+
+ return m_Ar;
+ }
+
+ PASTIS_INLINE
+ const NodeValue<const Rd>
+ computeBr(const NodeValuePerCell<Rdd>& Ajr,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const CellValue<const Rd>& uj,
+ const CellValue<const double>& pj) {
+ const auto& node_to_cell_matrix
+ = m_connectivity.nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells
+ = m_connectivity.nodeLocalNumbersInTheirCells();
+
+ parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+ Rd& br = m_br[r];
+ br = zero;
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells
+ = node_local_numbers_in_their_cells.itemValues(r);
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ br += Ajr(J,R)*uj[J] + pj[J]*Cjr(J,R);
+ }
+ });
+
+ return m_br;
+ }
+
+ void applyBoundaryConditions()
+ {
+ for (const auto& handler : m_boundary_condition_list) {
+ switch (handler.boundaryCondition().type()) {
+ case BoundaryCondition::normal_velocity: {
+ perr() << __FILE__ << ':' << __LINE__ << ": normal_velocity BC NIY\n";
+ std::exit(0);
+ break;
+ }
+ case BoundaryCondition::velocity: {
+ perr() << __FILE__ << ':' << __LINE__ << ": velocity BC NIY\n";
+ std::exit(0);
+ break;
+ }
+ case BoundaryCondition::pressure: {
+ // const PressureBoundaryCondition& pressure_bc
+ // = dynamic_cast<const PressureBoundaryCondition&>(handler.boundaryCondition());
+ perr() << __FILE__ << ':' << __LINE__ << ": pressure BC NIY\n";
+ std::exit(0);
+ break;
+ }
+ case BoundaryCondition::symmetry: {
+ const SymmetryBoundaryCondition<dimension>& symmetry_bc
+ = dynamic_cast<const SymmetryBoundaryCondition<dimension>&>(handler.boundaryCondition());
+ const Rd& n = symmetry_bc.outgoingNormal();
+
+ const Rdd I = identity;
+ const Rdd nxn = tensorProduct(n,n);
+ const Rdd P = I-nxn;
+
+ const Array<const NodeId>& node_list
+ = symmetry_bc.nodeList();
+ parallel_for(symmetry_bc.numberOfNodes(), PASTIS_LAMBDA(const int& r_number) {
+ const NodeId r = node_list[r_number];
+
+ m_Ar[r] = P*m_Ar[r]*P + nxn;
+ m_br[r] = P*m_br[r];
+ });
+ break;
+ }
+ }
+ }
+ }
+
+ NodeValue<Rd>
+ computeUr(const NodeValue<const Rdd>& Ar,
+ const NodeValue<const Rd>& br)
+ {
+ inverse(Ar, m_inv_Ar);
+ const NodeValue<const Rdd> invAr = m_inv_Ar;
+ parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+ m_ur[r]=invAr[r]*br[r];
+ });
+
+ return m_ur;
+ }
+
+ void
+ computeFjr(const NodeValuePerCell<Rdd>& Ajr,
+ const NodeValue<const Rd>& ur,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const CellValue<const Rd>& uj,
+ const CellValue<const double>& pj)
+ {
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().cellToNodeMatrix();
+
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ for (size_t r=0; r<cell_nodes.size(); ++r) {
+ m_Fjr(j,r) = Ajr(j,r)*(uj[j]-ur[cell_nodes[r]])+pj[j]*Cjr(j,r);
+ }
+ });
+ }
+
+ void inverse(const NodeValue<const Rdd>& A,
+ NodeValue<Rdd>& inv_A) const
+ {
+ parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+ inv_A[r] = ::inverse(A[r]);
+ });
+ }
+
+ PASTIS_INLINE
+ void computeExplicitFluxes(const NodeValue<const Rd>& xr,
+ const CellValue<const Rd>& xj,
+ const CellValue<const double>& rhoj,
+ const CellValue<const Rd>& uj,
+ const CellValue<const double>& pj,
+ const CellValue<const double>& cj,
+ const CellValue<const double>& Vj,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const NodeValuePerCell<const double>& ljr,
+ const NodeValuePerCell<const Rd>& njr) {
+ const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
+ computeAjr(rhocj, Cjr, ljr, njr);
+
+ NodeValuePerCell<const Rdd> Ajr = m_Ajr;
+ const NodeValue<const Rdd> Ar = computeAr(Ajr);
+ const NodeValue<const Rd> br = computeBr(m_Ajr, Cjr, uj, pj);
+
+ this->applyBoundaryConditions();
+
+ NodeValue<Rd>& ur = m_ur;
+ ur = computeUr(Ar, br);
+ computeFjr(m_Ajr, ur, Cjr, uj, pj);
+ }
+
+ NodeValue<Rd> m_br;
+ NodeValuePerCell<Rdd> m_Ajr;
+ NodeValue<Rdd> m_Ar;
+ NodeValue<Rdd> m_inv_Ar;
+ NodeValuePerCell<Rd> m_Fjr;
+ NodeValue<Rd> m_ur;
+ CellValue<double> m_rhocj;
+ CellValue<double> m_Vj_over_cj;
+
+ CellValue<Rd> m_xg;
+
+ public:
+ VoronoiSolver(MeshData& mesh_data,
+ UnknownsType& unknowns,
+ const std::vector<BoundaryConditionHandler>& bc_list)
+ : m_mesh_data(mesh_data),
+ m_mesh(mesh_data.mesh()),
+ m_connectivity(m_mesh.connectivity()),
+ m_boundary_condition_list(bc_list),
+ m_br(m_connectivity),
+ m_Ajr(m_connectivity),
+ m_Ar(m_connectivity),
+ m_inv_Ar(m_connectivity),
+ m_Fjr(m_connectivity),
+ m_ur(m_connectivity),
+ m_rhocj(m_connectivity),
+ m_Vj_over_cj(m_connectivity),
+ m_xg(m_connectivity)
+ {
+ ;
+ }
+
+ PASTIS_INLINE
+ double voronoi_dt(const CellValue<const double>& Vj,
+ const CellValue<const double>& cj) const
+ {
+ const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().cellToNodeMatrix();
+
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ double S = 0;
+ for (size_t r=0; r<cell_nodes.size(); ++r) {
+ S += ljr(j,r);
+ }
+ m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]);
+ });
+
+ return ReduceMin(m_Vj_over_cj);
+ }
+
+ void computeNextStep(const double& t, const double& dt,
+ UnknownsType& unknowns)
+ {
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<Rd>& uj = unknowns.uj();
+ CellValue<double>& Ej = unknowns.Ej();
+
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& cj = unknowns.cj();
+
+ const CellValue<const Rd>& xj = m_mesh_data.xj();
+
+ static bool init_xg=false;
+ if (not init_xg) {
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ m_xg[j] = xj[j];
+ });
+ init_xg = true;
+ }
+
+ const CellValue<const double>& Vj = m_mesh_data.Vj();
+ const NodeValuePerCell<const Rd>& Cjr = m_mesh_data.Cjr();
+ const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
+ const NodeValuePerCell<const Rd>& njr = m_mesh_data.njr();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
+
+ // computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr);
+
+ CellValue<Rd> uj_star(m_mesh.connectivity());
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ if ((j != 0) and ( j != m_mesh.numberOfCells()-1)) {
+ CellId jp1 = j+1;
+ CellId jm1 = j-1;
+ double rhoc_jp1 = rhoj[jp1]*cj[jp1];
+ double rhoc_jm1 = rhoj[jm1]*cj[jm1];
+ uj_star[j]
+ = (1./ (rhoc_jp1+rhoc_jm1)) * (rhoc_jp1 * uj[jp1] + rhoc_jm1*uj[jm1])
+ + (pj[jm1] - pj[jp1]) / (rhoc_jp1+rhoc_jm1);
+ } else {
+ uj_star[j] = zero;
+ }
+ });
+
+ // for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
+ // std::cout << "uj_star[" << j << "]= " << uj_star[j] << "\n";
+ // }
+
+
+ CellValue<double> pj_star(m_mesh.connectivity());
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ if ((j > 0) and ( j < m_mesh.numberOfCells()-1)) {
+ CellId jm1 = j-1;
+ double rhoc_jm1 = rhoj[jm1]*cj[jm1];
+ pj_star[j] = rhoc_jm1*(uj[jm1]-uj_star[j])[0] + pj[jm1];
+ } else {
+ pj_star[j] = pj[j];
+ }
+ });
+
+ const CellValue<const double> inv_mj = unknowns.invMj();
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ if ((j > 0) and ( j < m_mesh.numberOfCells()-1)) {
+ CellId jp1 = j+1;
+ CellId jm1 = j-1;
+ uj[j] -= (dt*inv_mj[j]) * 0.5*(pj_star[jp1]-pj_star[jm1]);
+ Ej[j] -= (dt*inv_mj[j]) * 0.5*(pj_star[jp1]*uj_star[jp1][0]-pj_star[jm1]*uj_star[jm1][0]);
+ }
+ });
+
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ m_xg[j] += dt*uj_star[j];
+ });
+
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+ ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
+ });
+
+ const auto& node_to_cell_matrix
+ = m_connectivity.nodeToCellMatrix();
+
+ NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
+ parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r){
+ const auto& node_to_cell = node_to_cell_matrix[r];
+
+ if (node_to_cell.size() == 2) {
+ Rd new_xr = zero;
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ new_xr += m_xg[J];
+ }
+ mutable_xr[r] = 0.5*new_xr;
+ }
+ });
+ m_mesh_data.updateAllData();
+
+ const CellValue<const double> mj = unknowns.mj();
+ parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+ rhoj[j] = mj[j]/Vj[j];
+ });
+
+
+ // for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
+ // std::cout << "rho[" << j << "]= " << rhoj[j] << "\t ";
+ // std::cout << "u[" << j << "]= " << uj[j] << "\t ";
+ // std::cout << "E[" << j << "]= " << Ej[j] << "\t ";
+ // std::cout << "e[" << j << "]= " << ej[j] << "\n";
+ // }
+ // std::cout << std::flush;
+ // std::cout << "enter value to continue: ";
+ // int i;
+ // std::cin >> i;
+ }
+};
+
+#endif // VORONOI_SOLVER_HPP