From 45d2aeb922bcd20ba7c76a85d4998536c0f9ad53 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Tue, 22 Dec 2020 10:11:25 +0100
Subject: [PATCH] Add an output frequency and cleaning
---
src/language/algorithms/ParabolicHeat.cpp | 22 +++++++++++-----
src/language/algorithms/ParabolicHeat.hpp | 3 ++-
src/language/modules/SchemeModule.cpp | 11 ++++----
src/scheme/ScalarDiamondScheme.hpp | 32 +++--------------------
4 files changed, 27 insertions(+), 41 deletions(-)
diff --git a/src/language/algorithms/ParabolicHeat.cpp b/src/language/algorithms/ParabolicHeat.cpp
index 4a0857bfa..e4527a682 100644
--- a/src/language/algorithms/ParabolicHeat.cpp
+++ b/src/language/algorithms/ParabolicHeat.cpp
@@ -32,7 +32,8 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
const double& Tf,
- const double& dt)
+ const double& dt,
+ const double& output_frequency)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
@@ -45,7 +46,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
return ItemType::node;
}
}();
- // HERE build CLs, time loop, main outputs
+
using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition,
SymmetryBoundaryCondition>;
@@ -244,7 +245,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
NodeValue<double> Tr(mesh->connectivity());
{
- VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01);
+ VTKWriter vtk_writer("Tanalytic_" + std::to_string(Dimension), 0.01);
vtk_writer.write(mesh,
{NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
@@ -256,12 +257,16 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
double time = 0;
Assert(dt > 0, "The time-step have to be positive!");
const CellValue<const double> primal_Vj = mesh_data.Vj();
+ double output_time = 0;
do {
{
- VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension) + std::to_string(time), 0.01);
+ if (time >= output_time) {
+ VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension) + "_time_" + std::to_string(time), 0.01);
- vtk_writer.write(mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj}}, 0,
- true); // forces last output
+ vtk_writer.write(mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj}}, 0,
+ true); // forces last output
+ output_time += output_frequency;
+ }
}
double deltat = std::min(dt, Tf - time);
std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
@@ -311,7 +316,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((face_error, face_error)) << "\n";
{
- VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
+ VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension) + "_final", 0.01);
vtk_writer.write(mesh,
{NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj},
@@ -449,6 +454,7 @@ template ParabolicHeatScheme<1>::ParabolicHeatScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
const double&,
+ const double&,
const double&);
template ParabolicHeatScheme<2>::ParabolicHeatScheme(
@@ -458,6 +464,7 @@ template ParabolicHeatScheme<2>::ParabolicHeatScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
const double&,
+ const double&,
const double&);
template ParabolicHeatScheme<3>::ParabolicHeatScheme(
@@ -467,4 +474,5 @@ template ParabolicHeatScheme<3>::ParabolicHeatScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
const double&,
+ const double&,
const double&);
diff --git a/src/language/algorithms/ParabolicHeat.hpp b/src/language/algorithms/ParabolicHeat.hpp
index 92df71a53..b5a2db881 100644
--- a/src/language/algorithms/ParabolicHeat.hpp
+++ b/src/language/algorithms/ParabolicHeat.hpp
@@ -25,7 +25,8 @@ class ParabolicHeatScheme
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
const double& final_time,
- const double& dt);
+ const double& dt,
+ const double&);
};
#endif // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ef8474bf3..e32f667e7 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -267,27 +267,28 @@ SchemeModule::SchemeModule()
void(std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
- const double&, const double&)>>(
+ const double&, const double&, const double&)>>(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
- const FunctionSymbolId& f_id, const double& final_time, const double& dt) -> void {
+ const FunctionSymbolId& f_id, const double& final_time, const double& dt,
+ const double& output_frequency) -> void {
switch (p_mesh->dimension()) {
case 1: {
ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id,
- f_id, final_time, dt};
+ f_id, final_time, dt, output_frequency};
break;
}
case 2: {
ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id,
- f_id, final_time, dt};
+ f_id, final_time, dt, output_frequency};
break;
}
case 3: {
ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id,
- f_id, final_time, dt};
+ f_id, final_time, dt, output_frequency};
break;
}
default: {
diff --git a/src/scheme/ScalarDiamondScheme.hpp b/src/scheme/ScalarDiamondScheme.hpp
index 87e3ae320..7259c8442 100644
--- a/src/scheme/ScalarDiamondScheme.hpp
+++ b/src/scheme/ScalarDiamondScheme.hpp
@@ -15,7 +15,6 @@
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItem.hpp>
-#include <output/VTKWriter.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
@@ -139,11 +138,6 @@ class ScalarDiamondScheme
~SymmetryBoundaryCondition() = default;
};
- // class DirichletBoundaryCondition;
- // class FourierBoundaryCondition;
- // class NeumannBoundaryCondition;
- // class SymmetryBoundaryCondition;
-
public:
ScalarDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
@@ -173,7 +167,6 @@ class ScalarDiamondScheme
BoundaryConditionList boundary_condition_list;
- std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
for (const auto& bc_descriptor : bc_descriptor_list) {
@@ -457,13 +450,6 @@ class ScalarDiamondScheme
diamond_mesh_data
.xj());
- VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
-
- vtk_writer.write(diamond_mesh,
- {NamedItemValue{"kappaj", dual_kappaj}, NamedItemValue{"coords", diamond_mesh->xr()},
- NamedItemValue{"Vj", diamond_mesh_data.Vj()}, NamedItemValue{"xj", diamond_mesh_data.xj()}},
- 0, true); // forces last output
-
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
const FaceValue<const double> mes_l = [&] {
@@ -574,27 +560,18 @@ class ScalarDiamondScheme
FaceValue<double> computed_alpha_l{mesh->connectivity()};
mapper->fromDualCell(alpha_j, computed_alpha_l);
- VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01);
-
- vtk_writer.write(diamond_mesh,
- {NamedItemValue{"alphaj", alpha_j}, NamedItemValue{"xj", diamond_mesh_data.xj()},
- NamedItemValue{"Vl", diamond_mesh_data.Vj()}, NamedItemValue{"|l|", dual_mes_l_j}},
- 0,
- true); // forces last output
-
return computed_alpha_l;
}();
double lambda = (Tf == 0) ? 0 : 1;
double time_factor = (lambda == 0) ? 1 : dt;
- std::cout << " lambda = " << lambda << " time_factor =" << time_factor << "\n";
+
const CellValue<const double> primal_Vj = mesh_data.Vj();
SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
- // const double beta_l = 1. / Dimension * alpha_l[face_id] * mes_l[face_id];
- const double beta_l = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
+ const double beta_l = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
const CellId cell_id1 = primal_face_to_cell[i_cell];
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
@@ -693,7 +670,7 @@ class ScalarDiamondScheme
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- b[face_dof_number[face_id]] += value_list[i_face]; // sign
+ b[face_dof_number[face_id]] += value_list[i_face];
}
}
},
@@ -711,8 +688,7 @@ class ScalarDiamondScheme
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
FaceId face_id = face_list[i_face];
- b[face_dof_number[face_id]] +=
- mes_l[face_id] * value_list[i_face]; // + lambda * Temperature_face[face_id]); // sign
+ b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];
}
}
},
--
GitLab