diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp index 7ac8e85f21026ff9c6f83b26d5b949d7590972e8..e6d79c1981300a661a21c499388e37884f0c32db 100644 --- a/src/scheme/FluxingAdvectionSolver.cpp +++ b/src/scheme/FluxingAdvectionSolver.cpp @@ -98,8 +98,10 @@ FluxingAdvectionSolver<2>::computeFluxVolume() const if (m_new_mesh->shared_connectivity() != m_old_mesh->shared_connectivity()) { throw NormalError("Old and new meshes must share the same connectivity"); } - MeshDataType& old_mesh_data = MeshDataManager::instance().getMeshData(*m_old_mesh); - MeshDataType& new_mesh_data = MeshDataManager::instance().getMeshData(*m_new_mesh); + std::cout << " CARRE " + << "\n"; + // MeshDataType& old_mesh_data = MeshDataManager::instance().getMeshData(*m_old_mesh); + // MeshDataType& new_mesh_data = MeshDataManager::instance().getMeshData(*m_new_mesh); const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix(); FaceValue<double> flux_volume(m_new_mesh->connectivity()); parallel_for( @@ -111,6 +113,8 @@ FluxingAdvectionSolver<2>::computeFluxVolume() const const Rd x3 = m_new_mesh->xr()[face_to_node[0]]; TinyMatrix<2> M(x2[0] - x0[0], x3[0] - x1[0], x2[1] - x0[1], x3[1] - x1[1]); flux_volume[face_id] = 0.5 * det(M); + // std::cout << " x1 " << x1 << " x0 " << x0 << " x3 " << x3 << " x2 " << x2 << " flux volume " + // << flux_volume[face_id] << "\n"; }); return flux_volume; } @@ -122,17 +126,84 @@ FluxingAdvectionSolver<3>::computeFluxVolume() const throw NotImplementedError("ViensViensViens"); } +template <typename MeshType> +auto +calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh, + [[maybe_unused]] const FaceValue<double>& fluxing_volumes) +{ + constexpr size_t Dimension = MeshType::Dimension; + const FaceValuePerCell<const bool> cell_face_is_reversed = old_mesh->connectivity().cellFaceIsReversed(); + const auto cell_to_face_matrix = old_mesh->connectivity().cellToFaceMatrix(); + const CellValue<double> total_negative_flux(old_mesh->connectivity()); + total_negative_flux.fill(0); + parallel_for( + old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const auto& cell_to_face = cell_to_face_matrix[cell_id]; + for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) { + FaceId face_id = cell_to_face[i_face]; + double flux = fluxing_volumes[face_id]; + if (cell_face_is_reversed(cell_id, i_face)) { + flux = -flux; + } + if (flux < 0) { + total_negative_flux[cell_id] += flux; + } + } + std::cout << " cell_id " << cell_id << " total_negative_flux " << total_negative_flux[cell_id] << "\n"; + }); + MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*old_mesh); + const CellValue<const double> Vj = mesh_data.Vj(); + const CellValue<size_t> ratio(old_mesh->connectivity()); + parallel_for( + old_mesh->numberOfCells(), + PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(abs(total_negative_flux[cell_id]) / Vj[cell_id]); }); + size_t number_of_cycle = max(ratio); + std::cout << " number_of_cycle " << number_of_cycle << "\n"; + return number_of_cycle; +} + template <typename MeshType, typename DataType> auto remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh, [[maybe_unused]] const FaceValue<double>& fluxing_volumes, + const size_t num, const DiscreteFunctionP0<MeshType::Dimension, const DataType>& old_q) { constexpr size_t Dimension = MeshType::Dimension; // const Connectivity<Dimension>& connectivity = new_mesh->connectivity(); - + const FaceValuePerCell<const bool> cell_face_is_reversed = new_mesh->connectivity().cellFaceIsReversed(); DiscreteFunctionP0<Dimension, DataType> new_q(new_mesh, copy(old_q.cellValues())); - + const auto cell_to_face_matrix = new_mesh->connectivity().cellToFaceMatrix(); + const auto face_to_cell_matrix = new_mesh->connectivity().faceToCellMatrix(); + for (size_t jstep = 0; jstep < num; ++jstep) { + std::cout << " step " << jstep << "\n"; + parallel_for( + new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const auto& cell_to_face = cell_to_face_matrix[cell_id]; + for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) { + FaceId face_id = cell_to_face[i_face]; + double flux = fluxing_volumes[face_id]; + if (cell_face_is_reversed(cell_id, i_face)) { + flux = -flux; + } + const auto& face_to_cell = face_to_cell_matrix[face_id]; + if (face_to_cell.size() == 1) { + continue; + } + CellId other_cell_id = face_to_cell[0]; + if (other_cell_id == cell_id) { + other_cell_id = face_to_cell[1]; + } + DataType fluxed_q = old_q[cell_id]; + if (flux > 0) { + fluxed_q = old_q[other_cell_id]; + } + fluxed_q *= flux / num; + new_q[cell_id] += fluxed_q; + } + // std::cout << " old q " << old_q[cell_id] << " new q " << new_q[cell_id] << "\n"; + }); + } return new_q; } @@ -140,6 +211,7 @@ template <typename MeshType> auto remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh, [[maybe_unused]] const FaceValue<double>& fluxing_volumes, + [[maybe_unused]] const size_t num, const DiscreteFunctionP0Vector<MeshType::Dimension, const double>& old_q) { constexpr size_t Dimension = MeshType::Dimension; @@ -244,21 +316,7 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh, switch (old_mesh->dimension()) { case 1: { - constexpr size_t Dimension = 1; - using MeshType = Mesh<Connectivity<Dimension>>; - - const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh); - - const std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh); - - // FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0); - - FaceValue<double> fluxing_volumes(old_mesh0->connectivity()); - fluxing_volumes.fill(0); - - return std::make_shared<DiscreteFunctionVariant>( - remapUsingFluxing(new_mesh0, fluxing_volumes, - remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>())); + throw NormalError("Not yet implemented in 1d"); } case 2: { constexpr size_t Dimension = 2; @@ -271,14 +329,10 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh, FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0); FaceValue<double> fluxing_volumes = solver.computeFluxVolume(); - //(new_mesh0->connectivity()); - // fluxing_volumes.fill(0); - + size_t number_of_cycles = calculateRemapCycles(old_mesh0, fluxing_volumes); return std::make_shared<DiscreteFunctionVariant>( - remapUsingFluxing(new_mesh0, fluxing_volumes, + remapUsingFluxing(new_mesh0, fluxing_volumes, number_of_cycles, remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>())); - - // throw NotImplementedError("Fluxing advection solver not implemented in dimension 2"); } case 3: { throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");