From a143de698384b1aa89457bba0ccc0722de419ecd Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Tue, 21 Feb 2023 17:29:02 +0100
Subject: [PATCH] fix bugs and allow vectors and tensors
---
src/scheme/FluxingAdvectionSolver.cpp | 81 ++++++++++++++++++++-------
1 file changed, 60 insertions(+), 21 deletions(-)
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 6d54194f1..b655ff74f 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -126,10 +126,10 @@ FluxingAdvectionSolver<3>::computeFluxVolume() const
throw NotImplementedError("ViensViensViens");
}
-template <typename MeshType>
+template <typename MeshType, typename DataType>
auto
calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
- [[maybe_unused]] const FaceValue<double>& fluxing_volumes)
+ [[maybe_unused]] const FaceValue<DataType>& fluxing_volumes)
{
constexpr size_t Dimension = MeshType::Dimension;
const FaceValuePerCell<const bool> cell_face_is_reversed = old_mesh->connectivity().cellFaceIsReversed();
@@ -142,7 +142,7 @@ calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
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)) {
+ if (!cell_face_is_reversed(cell_id, i_face)) {
flux = -flux;
}
if (flux < 0) {
@@ -164,8 +164,9 @@ calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
template <typename MeshType, typename DataType>
auto
-remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
- [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
+remapUsingFluxing(const std::shared_ptr<const MeshType>& old_mesh,
+ const std::shared_ptr<const MeshType>& new_mesh,
+ const FaceValue<double>& fluxing_volumes,
const size_t num,
const DiscreteFunctionP0<MeshType::Dimension, const DataType>& old_q)
{
@@ -173,17 +174,26 @@ remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
// 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();
+ DiscreteFunctionP0<Dimension, DataType> previous_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();
+ MeshData<Dimension>& old_mesh_data = MeshDataManager::instance().getMeshData(*old_mesh);
+ const CellValue<const double> oldVj = old_mesh_data.Vj();
+ const CellValue<double> Vjstep(new_mesh->connectivity());
+ parallel_for(
+ new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ Vjstep[cell_id] = oldVj[cell_id];
+ new_q[cell_id] *= oldVj[cell_id];
+ });
for (size_t jstep = 0; jstep < num; ++jstep) {
- std::cout << " step " << jstep << "\n";
+ // 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)) {
+ if (!cell_face_is_reversed(cell_id, i_face)) {
flux = -flux;
}
const auto& face_to_cell = face_to_cell_matrix[face_id];
@@ -194,30 +204,49 @@ remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
if (other_cell_id == cell_id) {
other_cell_id = face_to_cell[1];
}
- DataType fluxed_q = old_q[cell_id];
+ DataType fluxed_q = previous_q[cell_id];
if (flux > 0) {
- fluxed_q = old_q[other_cell_id];
+ fluxed_q = previous_q[other_cell_id];
}
+ Vjstep[cell_id] += flux / num;
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";
});
+ parallel_for(
+ new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ previous_q[cell_id] = 1 / Vjstep[cell_id] * new_q[cell_id];
+ // std::cout << " old q " << old_q[cell_id] << " new q " << previous_q[cell_id] << "\n";
+ // std::cout << " old vj " << oldVj[cell_id] << " new Vj " << Vjstep[cell_id] << "\n";
+ });
}
+
+ MeshData<Dimension>& new_mesh_data = MeshDataManager::instance().getMeshData(*new_mesh);
+ const CellValue<const double> newVj = new_mesh_data.Vj();
+ parallel_for(
+ new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] = 1 / newVj[cell_id] * new_q[cell_id]; });
+ // for (CellId cell_id = 0; cell_id < new_mesh->numberOfCells(); ++cell_id) {
+ // if (abs(newVj[cell_id] - Vjstep[cell_id]) > 1e-15) {
+ // std::cout << " cell " << cell_id << " newVj " << newVj[cell_id] << " Vjstep " << Vjstep[cell_id] << " diff "
+ // << abs(newVj[cell_id] - Vjstep[cell_id]) << "\n";
+ // }
+ // }
return new_q;
}
-template <typename MeshType>
+template <typename MeshType, typename DataType>
auto
-remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
+remapUsingFluxing([[maybe_unused]] const std::shared_ptr<const MeshType>& old_mesh,
+ 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)
+ const DiscreteFunctionP0Vector<MeshType::Dimension, const DataType>& old_q)
{
constexpr size_t Dimension = MeshType::Dimension;
// const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
- DiscreteFunctionP0Vector<Dimension, double> new_q(new_mesh, copy(old_q.cellArrays()));
+ DiscreteFunctionP0Vector<Dimension, DataType> new_q(new_mesh, copy(old_q.cellArrays()));
throw NotImplementedError("DiscreteFunctionP0Vector");
@@ -306,11 +335,11 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
std::shared_ptr<const DiscreteFunctionVariant>
FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
- const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables)
+ const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variable)
{
- const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({remapped_variables});
+ const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({remapped_variable});
- if (not checkDiscretizationType({remapped_variables}, DiscreteFunctionType::P0)) {
+ if (not checkDiscretizationType({remapped_variable}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
@@ -330,9 +359,19 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
FaceValue<double> fluxing_volumes = solver.computeFluxVolume();
size_t number_of_cycles = calculateRemapCycles(old_mesh0, fluxing_volumes);
- return std::make_shared<DiscreteFunctionVariant>(
- remapUsingFluxing(new_mesh0, fluxing_volumes, number_of_cycles,
- remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>()));
+
+ DiscreteFunctionVariant new_variable = std::visit(
+ [&](auto&& variable) -> DiscreteFunctionVariant {
+ using DiscreteFunctionT = std::decay_t<decltype(variable)>;
+ if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
+ return remapUsingFluxing(old_mesh0, new_mesh0, fluxing_volumes, number_of_cycles, variable);
+ } else {
+ throw UnexpectedError("incompatible mesh types");
+ }
+ },
+ remapped_variable->discreteFunction());
+
+ return std::make_shared<DiscreteFunctionVariant>(new_variable);
}
case 3: {
throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");
--
GitLab