Skip to content
Snippets Groups Projects
Commit 5407378f authored by Alexiane Plessier's avatar Alexiane Plessier
Browse files

Add Implicit traffic algorithm [ci skip]

parent 06cd12a1
No related branches found
No related tags found
No related merge requests found
...@@ -29,6 +29,7 @@ ...@@ -29,6 +29,7 @@
#include <scheme/NumberedBoundaryDescriptor.hpp> #include <scheme/NumberedBoundaryDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp> #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/UpwindExplicitTrafficSolver.hpp> #include <scheme/UpwindExplicitTrafficSolver.hpp>
#include <scheme/UpwindImplicitTrafficSolver.hpp>
#include <scheme/PerfectGas.hpp> #include <scheme/PerfectGas.hpp>
...@@ -150,6 +151,23 @@ SchemeModule::SchemeModule() ...@@ -150,6 +151,23 @@ SchemeModule::SchemeModule()
return UpwindExplicitTrafficSolverHandler{tau, bc_descriptor_list}.solver().apply(tau); return UpwindExplicitTrafficSolverHandler{tau, bc_descriptor_list}.solver().apply(tau);
})); }));
this->_addBuiltinFunction(
"implicit_traffic_flow_upwind",
std::make_shared<BuiltinFunctionEmbedder<
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<
const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const double& dt)>>(
[](const std::shared_ptr<const IDiscreteFunction>& tau,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const double& dt) -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>>
{
return UpwindImplicitTrafficSolverHandler{tau, bc_descriptor_list, dt}.solver().apply(dt, tau);
}));
this->_addBuiltinFunction("upwind_explicit_traffic", this->_addBuiltinFunction("upwind_explicit_traffic",
std::make_shared<BuiltinFunctionEmbedder< std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>, void(std::shared_ptr<const IMesh>,
......
...@@ -7,4 +7,5 @@ add_library( ...@@ -7,4 +7,5 @@ add_library(
DiscreteFunctionUtils.cpp DiscreteFunctionUtils.cpp
ImplicitAcousticSolver.cpp ImplicitAcousticSolver.cpp
UpwindExplicitTrafficSolver.cpp UpwindExplicitTrafficSolver.cpp
UpwindImplicitTrafficSolver.cpp
PerfectGas.cpp) PerfectGas.cpp)
...@@ -215,7 +215,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final ...@@ -215,7 +215,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
{ {
NodeValue<Rd> computed_ur(mesh.connectivity()); NodeValue<Rd> computed_ur(mesh.connectivity());
parallel_for( parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = -upwind_flux[r]; }); m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = -upwind_flux[r]; });
return computed_ur; return computed_ur;
} }
...@@ -256,7 +256,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final ...@@ -256,7 +256,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
NodeValue<Rd> new_xr = copy(mesh->xr()); NodeValue<Rd> new_xr = copy(mesh->xr());
parallel_for( parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; }); m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr); std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
...@@ -267,7 +267,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final ...@@ -267,7 +267,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
CellValue<double> new_tau = copy(tau.cellValues()); CellValue<double> new_tau = copy(tau.cellValues());
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j]; const auto& cell_nodes = cell_to_node_matrix[j];
double flux_sum = 0; double flux_sum = 0;
for (size_t r = 0; r < cell_nodes.size(); ++r) { for (size_t r = 0; r < cell_nodes.size(); ++r) {
...@@ -373,16 +373,16 @@ UpwindExplicitTrafficSolverHandler::UpwindExplicitTrafficSolverHandler( ...@@ -373,16 +373,16 @@ UpwindExplicitTrafficSolverHandler::UpwindExplicitTrafficSolverHandler(
std::make_unique<UpwindExplicitTrafficSolver<1>>(i_mesh, tau, bc_descriptor_list); std::make_unique<UpwindExplicitTrafficSolver<1>>(i_mesh, tau, bc_descriptor_list);
break; break;
} }
case 2: { // case 2: {
m_upwind_explicit_traffic_solver = // m_upwind_explicit_traffic_solver =
std::make_unique<UpwindExplicitTrafficSolver<2>>(i_mesh, tau, bc_descriptor_list); // std::make_unique<UpwindExplicitTrafficSolver<2>>(i_mesh, tau, bc_descriptor_list);
break; // break;
} // }
case 3: { // case 3: {
m_upwind_explicit_traffic_solver = // m_upwind_explicit_traffic_solver =
std::make_unique<UpwindExplicitTrafficSolver<3>>(i_mesh, tau, bc_descriptor_list); // std::make_unique<UpwindExplicitTrafficSolver<3>>(i_mesh, tau, bc_descriptor_list);
break; // break;
} // }
default: { default: {
throw UnexpectedError("invalid mesh dimension"); throw UnexpectedError("invalid mesh dimension");
} }
......
#include <scheme/UpwindImplicitTrafficSolver.hpp>
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/FreeBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunction.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <variant>
#include <vector>
template <size_t Dimension>
class UpwindImplicitTrafficSolverHandler::UpwindImplicitTrafficSolver final
: public UpwindImplicitTrafficSolverHandler::IUpwindImplicitTrafficSolver
{
using Rdxd = TinyMatrix<Dimension>;
using Rd = TinyVector<Dimension>;
using MeshType = Mesh<Connectivity<Dimension>>;
using MeshDataType = MeshData<Dimension>;
using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>;
class FreeBoundaryCondition;
class VelocityBoundaryCondition;
using BoundaryCondition = std::variant<FreeBoundaryCondition, VelocityBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
BoundaryConditionList m_boundary_condition_list;
const MeshType& m_mesh;
NodeValue<const Rd> m_ur;
NodeValue<const Rd> m_upwind_flux;
CellValue<const Rd> m_Mj;
CellValue<const double>
_getRho(const DiscreteScalarFunction& tau)
{
CellValue<double> rho{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { rho[j] = 1. / tau[j]; });
return rho;
}
CellValue<const double>
_getU(const DiscreteScalarFunction& tau)
{
CellValue<double> u{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { u[j] = 1 - 1. / tau[j]; });
return u;
}
const CellValue<const double>
_getMj(const DiscreteScalarFunction& tau)
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> Mj{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { Mj[j] = Vj[j] * 1. / tau[j]; });
return Mj;
}
const CellValue<const double>
_getInv_Mj(const DiscreteScalarFunction& tau)
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> inv_Mj{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { inv_Mj[j] = tau[j] * 1. / Vj[j]; });
return inv_Mj;
}
// function f
Vector<double>
_getF(const Vector<double>& taun, const Vector<double>& tauk, const double dt)
{
Vector<double> function_f{m_mesh.numberOfCells()};
const auto& cell_to_node_matrix = m_mesh.cellToNodeMatrix();
const auto& node_to_cell_matrix = m_mesh.nodeToCellMatrix();
const CellValue<const double> inv_Mj = [&](const DiscreteScalarFunction& tau) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> computed_inv_Mj{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { inv_Mj[j] = tau[j] * 1. / Vj[j]; });
return computed_inv_Mj;
}();
CellValue<double> computed_f(m_mesh.connectivity());
for (CellId j = 0; j < m_mesh.numberOfCells() - 1; ++j) {
const auto& cell_node = cell_to_node_matrix[j];
NodeId r = cell_node[1];
const auto& node_cell = node_to_cell_matrix[r];
CellId j1 = node_cell[1];
computed_f[j] = tauk[j] - taun[j] + (dt * inv_Mj[j]) / tauk[j1] - (dt * inv_Mj[j]) / tauk[j];
}
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId& node_id = node_list[i_node];
const auto& node_cell = node_to_cell_matrix[node_id];
CellId j = node_cell[0];
computed_f[j] =
tauk[j] - taun[j] + (dt * inv_Mj[j]) * (1 - value_list[i_node][0]) - (dt * inv_Mj[j]) / tauk[j];
}
}
},
boundary_condition);
}
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
size_t i = j;
function_f[i] = computed_f[j];
});
return function_f;
}
CRSMatrix<double>
_getGradientF(const Vector<double>& tauk, const double dt)
{
SparseMatrixDescriptor grad_f{m_mesh.numberOfCells()};
const auto& cell_to_node_matrix = m_mesh.cellToNodeMatrix();
const auto& node_to_cell_matrix = m_mesh.nodeToCellMatrix();
const CellValue<const double> inv_Mj = [&](const DiscreteScalarFunction& tau) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> computed_inv_Mj{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { inv_Mj[j] = tau[j] * 1. / Vj[j]; });
return computed_inv_Mj;
}();
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId& node_id = node_list[i_node];
const auto& node_cell = node_to_cell_matrix[node_id];
CellId j = node_cell[0];
size_t number_of_the_cell = j;
int row = number_of_the_cell;
int k = 0;
grad_f(k, row) -= (dt * inv_Mj[j]) / (tauk[j] * tauk[j]);
}
}
},
boundary_condition);
}
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
size_t number_of_the_cell = j;
int row = number_of_the_cell;
grad_f(row, row) = 1 + (dt * inv_Mj[j]) / (tauk[j] * tauk[j]);
}
for (CellId j = 0; j < m_mesh.numberOfCells() - 1; ++j) {
const auto& cell_node = cell_to_node_matrix[j];
NodeId r = cell_node[1];
const auto& node_cell = node_to_cell_matrix[r];
CellId j1 = node_cell[1];
size_t number_of_the_cell = j;
int row = number_of_the_cell;
grad_f(row, row + 1) = -dt * inv_Mj[j] / (tauk[j1] * tauk[j1]);
}
}
UpwindImplicitTrafficSolver(
const std::shared_ptr<const MeshType>& p_mesh,
const DiscreteScalarFunction& tau,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt)
: m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)},
m_mesh{*p_mesh}, // m_tau{0} pour vecteur initial de solution?
{
CellValue<double> taun = copy(tau); // on l'a deja dans las variables, mais peut etre faudrait le mettre en
// vecteur pour avoir les valeurs absolues, donc faire un gettau
CellValue<double> tauk = copy(taun);
int nb_iter = 0;
double norm_inf_sol;
CellValue<double> abs_taun = [&]() {
CellValue<double> computed_abs_taun(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { computed_abs_taun[j] = std::abs(taun[j]); });
return computed_abs_taun;
}();
double norm_inf_taun = max(abs_taun);
do {
std::cout << "iteration=" << nb_iter << '\n';
nb_iter++;
Vector<double> f = this->_getF(taun, tauk, dt);
CRSMatrix<double> gradient_f = this->_getGradientF(tauk, dt);
Vector<double> sol{m_mesh.numberOfCells()};
sol = 1;
LinearSolver solver;
solver.solveLocalSystem(gradient_f, sol, f);
Vector<double> tau_next = tauk - sol;
Vector<const double> abs_sol = [&]() {
Vector<double> compute_abs_sol{sol.size()};
parallel_for(
sol.size(), PUGS_LAMBDA(size_t i) { compute_abs_sol[i] = std::abs(sol[i]); });
return compute_abs_sol;
}();
norm_inf_sol = max(abs_sol);
tauk = tau_next;
std::cout << "ratio" << norm_inf_sol / norm_inf_taun << "\n";
} while ((norm_inf_sol > 1e-14 * norm_inf_taun) and (nb_iter < 10000));
m_tau = tauk;
}
public:
std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
apply(const double& dt,
const std::shared_ptr<const MeshType>& mesh,
const DiscreteFunctionP0<Dimension, double>& tau) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const NodeValue<Rd> new_position = [&]() {
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
NodeValue<Rd> x_next(mesh->connectivity());
NodeId r = 0;
const auto& node_cells = node_to_cell_matrix[r];
CellId j = node_cells[0];
x_next[r] = mesh->xr()[r][0] + dt * (1 - 1. / tau[j]);
for (NodeId r = 1; r < mesh->connectivity().numberOfNodes() - 1; ++r) {
const auto& node_cells = node_to_cell_matrix[r];
CellId j = node_cells[1];
x_next[r] = mesh->xr()[r][0] + dt * (1 - 1. / tau[j]);
}
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId& node_id = node_list[i_node];
x_next[node_id] = mesh->xr()[node_id][0] + dt * value_list[i_node][0];
}
}
},
boundary_condition);
}
return x_next;
}();
NodeValue<Rd> new_xr = copy(mesh->xr());
parallel_for(
mesh->connectivity().numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] = new_position[r]; });
mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
for (CellId j = 0; j < mesh->connectivity().numberOfCells(); ++j) {
tau[j] = tau_next[j];
// std::cout << '\n' << "tauj[" << j << "]=" << tauj[j] << '\n';
// std::cout << '\n' << "rhoj[" << j << "]=" << rhoj[j] << '\n';
// std::cout << '\n' << "uj[" << j << "]=" << uj[j] << '\n';
}
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
CellValue<double> new_tau = copy(tau.cellValues());
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_tau)};
}
std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>>
apply(const double& dt, const std::shared_ptr<const IDiscreteFunction>& tau) const
{
std::shared_ptr i_mesh = getCommonMesh({tau});
if (not i_mesh) {
throw NormalError("discrete function is not defined on the same mesh");
}
if (not checkDiscretizationType({tau}, DiscreteFunctionType::P0)) {
throw NormalError("traffic solver expects P0 functions");
}
return this->apply(dt, std::dynamic_pointer_cast<const MeshType>(i_mesh),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(tau));
}
UpwindImplicitTrafficSolver(
const std::shared_ptr<const IMesh>& mesh,
const std::shared_ptr<const IDiscreteFunction>& tau,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt)
: UpwindImplicitTrafficSolver(std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(tau),
bc_descriptor_list,
dt)
{}
UpwindImplicitTrafficSolver() = default;
UpwindImplicitTrafficSolver(UpwindImplicitTrafficSolver&&) = default;
~UpwindImplicitTrafficSolver() = default;
};
template <size_t Dimension>
class UpwindImplicitTrafficSolverHandler::UpwindImplicitTrafficSolver<Dimension>::VelocityBoundaryCondition
{
private:
const Array<const TinyVector<Dimension>> m_value_list;
const Array<const NodeId> m_node_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_node_list;
}
const Array<const TinyVector<Dimension>>&
valueList() const
{
return m_value_list;
}
VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
: m_value_list{value_list}, m_node_list{node_list}
{}
~VelocityBoundaryCondition() = default;
};
template <size_t Dimension>
class UpwindImplicitTrafficSolverHandler::UpwindImplicitTrafficSolver<Dimension>::FreeBoundaryCondition
{
private:
const Array<const NodeId> m_node_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_node_list;
}
FreeBoundaryCondition(const Array<const NodeId>& node_list) : m_node_list{node_list} {}
~FreeBoundaryCondition() = default;
};
UpwindImplicitTrafficSolverHandler::UpwindImplicitTrafficSolverHandler(
const std::shared_ptr<const IDiscreteFunction>& tau,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt)
{
std::shared_ptr i_mesh = getCommonMesh({tau});
if (not i_mesh) {
throw NormalError("discrete function is not defined on the same mesh");
}
if (not checkDiscretizationType({tau}, DiscreteFunctionType::P0)) {
throw NormalError("traffic solver expects P0 functions");
}
switch (i_mesh->dimension()) {
case 1: {
m_upwind_implicit_traffic_solver =
std::make_unique<UpwindImplicitTrafficSolver<1>>(i_mesh, tau, bc_descriptor_list, dt);
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
#ifndef UPWIND_IMPLICIT_TRAFFIC_SOLVER_HPP
#define UPWIND_IMPLICIT_TRAFFIC_SOLVER_HPP
#include <memory>
#include <tuple>
#include <vector>
class IDiscreteFunction;
class IBoundaryConditionDescriptor;
class IMesh;
class UpwindImplicitTrafficSolverHandler
{
private:
struct IUpwindImplicitTrafficSolver
{
virtual std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>> apply(
const double& dt,
const std::shared_ptr<const IDiscreteFunction>& tau) const = 0;
IUpwindImplicitTrafficSolver() = default;
IUpwindImplicitTrafficSolver(IUpwindImplicitTrafficSolver&&) = default;
IUpwindImplicitTrafficSolver& operator=(IUpwindImplicitTrafficSolver&&) = default;
virtual ~IUpwindImplicitTrafficSolver() = default;
};
template <size_t Dimension>
class UpwindImplicitTrafficSolver;
std::unique_ptr<IUpwindImplicitTrafficSolver> m_upwind_implicit_traffic_solver;
public:
const IUpwindImplicitTrafficSolver&
solver() const
{
return *m_upwind_implicit_traffic_solver;
}
UpwindImplicitTrafficSolverHandler(
const std::shared_ptr<const IDiscreteFunction>& tau,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt);
};
#endif // UPWIND_IMPLICIT_TRAFFIC_SOLVER_HPP
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment