Skip to content
Snippets Groups Projects
Commit f2c3a45b authored by Axelle Drouard's avatar Axelle Drouard
Browse files

Change arguments of "nodalheat" for implicit scheme.

parent fc63eeaa
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/language/modules/SchemeModule.cpp
+ 18
− 16
View file @ f2c3a45b
......@@ -603,20 +603,22 @@ SchemeModule::SchemeModule()
}
));
this->_addBuiltinFunction("nodalheat",
std::make_shared<BuiltinFunctionEmbedder<
std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
this
->_addBuiltinFunction("nodalheat",
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
const std::vector<std::shared_ptr<
const IBoundaryConditionDescriptor>>&)>>(
const std::vector<
std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const std::shared_ptr<const IDiscreteFunction>&, const double&)>>(
[](const std::shared_ptr<const IDiscreteFunction>& mub_dual,
const std::shared_ptr<const IDiscreteFunction>& mu_dual,
[](const std::shared_ptr<const IDiscreteFunction>& kappa,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
return ScalarNodalSchemeHandler{mub_dual, mu_dual, f, bc_descriptor_list}.solution();
bc_descriptor_list,
const std::shared_ptr<const IDiscreteFunction>& P,
const double& dt) -> const std::shared_ptr<const IDiscreteFunction> {
return ScalarNodalSchemeHandler{kappa, f, bc_descriptor_list, P, dt}.solution();
}
));
......
src/scheme/ScalarNodalScheme.cpp
+ 40
− 79
View file @ f2c3a45b
......@@ -143,10 +143,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_b,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& P,
const double& dt)
{
using BoundaryCondition =
std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>;
......@@ -409,7 +410,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
{
CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_b->cellValues();
const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
......@@ -429,24 +429,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return kappa;
}();
const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
kappa[node_id] = zero;
const auto& node_to_cell = node_to_cell_matrix[node_id];
double weight = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
double local_weight = 1. / l2Norm(xr[node_id] - xj[cell_id]);
kappa[node_id] += local_weight * cell_kappajb[cell_id];
weight += local_weight;
}
kappa[node_id] = 1. / weight * kappa[node_id];
});
return kappa;
}();
const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
parallel_for(
......@@ -474,17 +456,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return kappa_invBeta;
}();
const NodeValue<const TinyMatrix<Dimension>> kapparb_invBetar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (not node_is_corner[node_id]) {
kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]);
}
});
return kappa_invBeta;
}();
const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] {
NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
parallel_for(
......@@ -538,17 +509,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return kappa_invBeta;
}();
const NodeValue<const TinyMatrix<Dimension>> corner_kapparb_invBetar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (node_is_corner[node_id]) {
kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(corner_betar[node_id]);
}
});
return kappa_invBeta;
}();
const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
parallel_for(
......@@ -570,8 +530,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (is_boundary_node[node_id]) {
compute_v[node_id] = 1. / l2Norm(node_kapparb[node_id] * exterior_normal[node_id]) *
node_kapparb[node_id] * exterior_normal[node_id];
compute_v[node_id] = 1. / l2Norm(node_kappar[node_id] * exterior_normal[node_id]) * node_kappar[node_id] *
exterior_normal[node_id];
}
});
return compute_v;
......@@ -634,7 +594,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const TinyMatrix<Dimension> I = identity;
const TinyVector<Dimension> n = exterior_normal[node_id];
const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
const TinyMatrix<Dimension> P = I - nxn;
const TinyMatrix<Dimension> Q = I - nxn;
for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
const CellId cell_id_j = node_to_cell[i_cell_j];
......@@ -647,7 +607,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
S(cell_id_j, cell_id_k) +=
dot(kappar_invBetar[node_id] *
(Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
P * Cjr(cell_id_k, i_node_j));
Q * Cjr(cell_id_k, i_node_j));
}
}
} else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) {
......@@ -679,18 +639,25 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}
for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) {
for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) {
S(cell_id_j, cell_id_k) *= dt / Vj[cell_id_j];
}
};
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
S(cell_id, cell_id) += Vj[cell_id];
S(cell_id, cell_id) += 1;
}
CellValue<const double> fj = f->cellValues();
CellValue<const double> Ph = P->cellValues();
CRSMatrix A{S.getCRSMatrix()};
Vector<double> b{mesh->numberOfCells()};
b = zero;
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
b[cell_id] = fj[cell_id] * Vj[cell_id];
b[cell_id] = fj[cell_id];
};
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
......@@ -699,16 +666,16 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const TinyMatrix<Dimension> I = identity;
const TinyVector<Dimension> n = exterior_normal[node_id];
const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
const TinyMatrix<Dimension> P = I - nxn;
const TinyMatrix<Dimension> Q = I - nxn;
if ((node_is_neumann[node_id]) && (not node_is_dirichlet[node_id])) {
if ((is_boundary_node[node_id]) and (not node_is_corner[node_id])) {
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node = node_local_number_in_its_cells[i_cell];
b[cell_id] += node_boundary_values[node_id] *
(1. / (sum_theta[node_id] * l2Norm(node_kapparb[node_id] * n)) *
dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) +
b[cell_id] += 1. / Vj[cell_id] * node_boundary_values[node_id] *
(1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
dot(Cjr(cell_id, i_node), n));
}
} else if (node_is_corner[node_id]) {
......@@ -719,7 +686,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
FaceId face_id = node_to_face[i_face];
sum_mes_l += mes_l[face_id];
}
b[cell_id] += 0.5 * node_boundary_values[node_id] * sum_mes_l;
b[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_values[node_id] * sum_mes_l;
}
} else if (node_is_dirichlet[node_id]) {
if (not node_is_corner[node_id]) {
......@@ -731,8 +698,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
b[cell_id_j] +=
dot(node_boundary_values[node_id] * kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
b[cell_id_j] += 1. / Vj[cell_id_j] *
dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
Cjr(cell_id_j, i_node_j));
}
}
......@@ -747,7 +714,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
b[cell_id_j] +=
dot(node_boundary_values[node_id] * corner_kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
1. / Vj[cell_id_j] *
dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
Cjr(cell_id_j, i_node_j));
}
}
......@@ -757,8 +725,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
Vector<double> T{mesh->numberOfCells()};
T = zero;
Vector<double> B{mesh->numberOfCells()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; });
LinearSolver solver;
solver.solveLocalSystem(A, T, b);
solver.solveLocalSystem(A, T, B);
m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
auto& solution = *m_solution;
......@@ -776,17 +748,18 @@ ScalarNodalSchemeHandler::solution() const
}
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& cell_k_b,
const std::shared_ptr<const IDiscreteFunction>& cell_k,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::shared_ptr<const IDiscreteFunction>& P,
const double& dt)
{
const std::shared_ptr i_mesh = getCommonMesh({cell_k, f});
if (not i_mesh) {
throw NormalError("primal discrete functions are not defined on the same mesh");
}
checkDiscretizationType({cell_k_b, cell_k, f}, DiscreteFunctionType::P0);
checkDiscretizationType({cell_k, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
......@@ -801,26 +774,14 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh,
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
bc_descriptor_list,
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(P), dt);
break;
}
case 3: {
using MeshType = Mesh<Connectivity<3>>;
using DiscreteTensorFunctionType = DiscreteFunctionP0<3, TinyMatrix<3>>;
using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>;
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
m_scheme =
std::make_unique<ScalarNodalScheme<3>>(mesh,
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
throw NormalError("The scheme is not implemented in dimension 3");
break;
}
default: {
......
src/scheme/ScalarNodalScheme.hpp
+ 4
− 3
View file @ f2c3a45b
......@@ -40,10 +40,11 @@ class ScalarNodalSchemeHandler
std::shared_ptr<const IDiscreteFunction> solution() const;
ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_k_bound,
const std::shared_ptr<const IDiscreteFunction>& cell_k,
ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_k,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::shared_ptr<const IDiscreteFunction>& P,
const double& dt);
~ScalarNodalSchemeHandler();
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment