1#ifndef DUNE_COPASI_MODEL_LOCAL_EQUATIONS_LOCAL_EQUATIONS_HH
2#define DUNE_COPASI_MODEL_LOCAL_EQUATIONS_LOCAL_EQUATIONS_HH
13#include <dune/pdelab/concepts/basis.hh>
14#include <dune/pdelab/common/container_entry.hh>
16#include <dune/geometry/dimension.hh>
18#include <dune/common/fvector.hh>
19#include <dune/common/indices.hh>
20#include <dune/common/overloadset.hh>
21#include <dune/common/parametertree.hh>
22#include <dune/common/tuplevector.hh>
24#include <function2/function2.hpp>
25#include <spdlog/spdlog.h>
37template<std::
size_t dim>
40 using Scalar = FieldVector<double, 1>;
41 using Vector = FieldVector<double, dim>;
42 using Tensor = FieldMatrix<double, dim, dim>;
44 using CompartmentPath = TypeTree::HybridTreePath<index_constant<0>, std::size_t, std::size_t>;
45 using MembranePath = TypeTree::HybridTreePath<index_constant<1>, std::size_t, std::size_t>;
62 assert(path[Indices::_1] == 0);
63 return lbasis.child(path[Indices::_2]);
70 return lbasis.child(path[Indices::_1]).child(path[Indices::_2]);
77 assert(path[Indices::_1] == 0);
78 return lbasis.child(path[Indices::_2]);
81 static const Concept::MembraneSubEntitiesLocalBasisNode
auto& path_to_local_basis_node(
83 const Concept::MultiCompartmentLocalBasisNode
auto& lbasis)
85 return lbasis.child(path[Indices::_1]).child(path[Indices::_2]);
88 enum class FactoryFalgs
98 struct CompartmentNode;
101 template<
class Signature>
106 : fu2::unique_function<Signature>{ std::move(callable) }
113 template<
class Signature>
118 : fu2::unique_function<Signature>{ std::move(callable) }
125 template<
class Signature>
131 template<
class Signature>
163 fu2::unique_function<Vector(Vector)
const noexcept>&& callable,
190 const PDELab::Concept::LocalBasis
auto& lbasis)
const
192 return path_to_local_basis_node(
path, lbasis.tree());
197 std::cout << fmt::format(
"Name: {}\n",
name);
198 std::cout <<
"\tPath: " <<
path << std::endl;
199 std::cout << fmt::format(
"\tValue: {}\n",
value[0]);
200 std::cout <<
"\tGradient: " <<
gradient << std::endl;
202 std::cout << fmt::format(
"\tReaction: {}\n",
reaction()[0]);
204 std::cout << fmt::format(
"\t\tJacobian wrt '{}': {}\n", jac.wrt.name, jac()[0]);
207 std::cout << fmt::format(
"\tStorage: {}\n",
storage()[0]);
209 std::cout << fmt::format(
"\tCross Diffusion wrt '{}': ", diffusion.wrt.name)
210 << diffusion(diffusion.wrt.gradient) << std::endl;
211 for (
const auto& jac : diffusion.compartment_jacobian)
212 std::cout << fmt::format(
"\t\tJacobian wrt '{}': ", jac.wrt.name) << jac(jac.wrt.gradient)
235 const PDELab::Concept::LocalBasis
auto& lbasis)
const
237 return path_to_local_basis_node(
path, lbasis.tree());
242 std::cout << fmt::format(
"Name: {}\n",
name);
243 std::cout <<
"\tPath: " <<
path << std::endl;
244 std::cout << fmt::format(
"\tValue: {}\n",
value[0]);
245 std::cout <<
"\tGradient: " <<
gradient << std::endl;
247 std::cout << fmt::format(
"\tReaction: {}\n",
reaction()[0]);
249 std::cout << fmt::format(
"\t\tJacobian wrt '{}': {}\n", jac.wrt.name, jac()[0]);
252 std::cout << fmt::format(
"\tStorage: {}\n",
storage()[0]);
254 std::cout << fmt::format(
"\tCross Diffusion wrt '{}': ", diffusion.wrt.name)
255 << diffusion(diffusion.wrt.gradient) << std::endl;
256 for (
const auto& jac : diffusion.compartment_jacobian)
257 std::cout << fmt::format(
"\t\tJacobian wrt '{}': ", jac.wrt.name) << jac(jac.wrt.gradient)
263 template<PDELab::Concept::LocalBasis LocalBasis,
class CellDataGr
idView =
typename LocalBasis::GlobalBasis::EntitySet,
class CellDataType =
double>
264 [[nodiscard]]
static std::unique_ptr<LocalEquations>
make(
265 const LocalBasis& lbasis,
266 const ParameterTree& eqs_cfg = {},
267 const ParameterTree& domain_cfg = {},
268 std::shared_ptr<const FunctorFactory<dim>> functor_factory =
nullptr,
269 std::shared_ptr<const CellData<CellDataGridView, CellDataType>> grid_cell_data =
nullptr,
272 TRACE_EVENT(
"dune",
"LocalEquations::make");
273 auto local_values = std::unique_ptr<LocalEquations>(
new LocalEquations{});
276 std::size_t compartment_count = 0;
277 std::size_t membrane_count = 0;
278 PDELab::forEachLeafNode(
280 overload([&](
const Concept::CompartmentScalarLocalBasisNode
auto&) { ++compartment_count; },
281 [&](
const Concept::MembraneScalarLocalBasisNode
auto&) { ++membrane_count; }));
282 local_values->_values.reserve(compartment_count + membrane_count);
283 local_values->_gradients.reserve(compartment_count + membrane_count);
285 local_values->initialize_nodes(
286 lbasis.tree(), [&](
auto path) { return lbasis.globalBasis().subSpace(path).name(); });
289 std::set<std::string_view> names;
290 PDELab::forEach(local_values->_nodes, [&](
auto& compartments_fncs) {
291 for (const auto& compartment_fncs : compartments_fncs)
292 for (const auto& component_fncs : compartment_fncs)
293 if (auto [it, inserted] = names.insert(component_fncs.name); not inserted)
294 throw format_exception(
295 InvalidStateException{},
"\tVariable with name '{}' is repeated", *it);
299 if (grid_cell_data) {
300 local_values->cell_values.resize( grid_cell_data->size() );
301 local_values->cell_mask.resize( grid_cell_data->size() );
302 local_values->cell_keys = grid_cell_data->keys();
305 if ((opts.
any() or domain_cfg.hasSub(
"deformation")) and not functor_factory)
307 "Equations cannot be configured without a functor factory");
310 local_values->configure_eqs(eqs_cfg, functor_factory, opts);
312 if (domain_cfg.hasSub(
"deformation"))
313 set_differentiable_function(*local_values,
315 local_values->domain_deformation,
316 "domain.deformation",
317 domain_cfg.sub(
"deformation"),
322 template<
class CellDataGr
idView,
class CellDataType =
double>
324 const PDELab::Concept::LocalBasis
auto& lbasis,
325 const ParameterTree& eqs_cfg,
326 const ParameterTree& domain_cfg,
336 FactoryFalgs::Reaction | FactoryFalgs::Diffusion | FactoryFalgs::Velocity |
337 FactoryFalgs::Outflow);
340 template<
class CellDataGr
idView,
class CellDataType =
double>
342 const PDELab::Concept::LocalBasis
auto& lbasis,
343 const ParameterTree& eqs_cfg,
344 const ParameterTree& domain_cfg,
354 FactoryFalgs::Storage);
362 return PDELab::containerEntry(_nodes, make_path(tree));
370 return PDELab::containerEntry(_nodes, make_path(tree)).value;
378 return PDELab::containerEntry(_nodes, make_path(tree)).gradient;
382 const std::function<
void(std::string_view,
383 const FieldVector<double, 1>&,
384 const FieldVector<double, dim>&)>& apply)
const override
386 PDELab::forEach(nodes(), [&](
auto& compartments) {
387 for (
auto& compartment_fncs : compartments)
388 for (
auto& component_fncs : compartment_fncs)
389 apply(component_fncs.name, component_fncs.value, component_fncs.gradient);
394 const std::function<
void(std::string_view,
395 const FieldVector<double, 1>&,
396 const FieldVector<double, dim - 1>&)>&)
const override
403 std::fill(begin(_values), end(_values), 0.);
404 std::fill(begin(_gradients), end(_gradients), 0.);
409 PDELab::forEach(_nodes, [&](
auto& compartments) {
410 for (
auto& compartment_fncs : compartments)
411 for (
auto& func : compartment_fncs)
416 const auto&
nodes()
const {
return _nodes; }
424 auto make_path(
const Tree& tree)
const
430 auto compartment_type = compartment_index(tree);
434 auto front_index = front(tree.orderingViewPath());
437 if constexpr (IsCompileTimeConstant<
decltype(front_index)>{}) {
438 assert(front_index == compartment_type);
439 return pop_front(tree.orderingViewPath());
442 return tree.orderingViewPath();
445 assert(mi.size() <= 2);
448 auto field_id = back(mi);
451 auto compartment_id = (mi.size() == 2) ? front(mi) : 0;
452 return TypeTree::treePath(compartment_type, compartment_id, field_id);
456 TupleVector<std::vector<std::vector<CompartmentNode>>, std::vector<std::vector<MembraneNode>>>
458 std::array<std::vector<std::string>, 2> _compartment_names;
460 std::vector<Scalar> _values;
461 std::vector<Vector> _gradients;
471 requires Concept::CompartmentLocalBasisNode<Tree> || Concept::MembraneLocalBasisNode<Tree>
472 void initialize_nodes(
const Tree& tree,
auto fname, std::size_t c = 0)
474 auto ct = compartment_index(tree);
475 _nodes[ct].resize(c + 1);
476 _nodes[ct][c].clear();
477 _compartment_names[ct].resize(c + 1);
478 for (std::size_t i = 0; i != tree.degree(); ++i) {
479 _nodes[ct][c].emplace_back(_values.emplace_back(),
480 _gradients.emplace_back(),
481 TypeTree::treePath(ct, c, i),
482 fname(tree.child(i).orderingViewPath()));
483 _compartment_names[ct][c] = fname(pop_back(tree.child(i).orderingViewPath()));
488 requires Concept::MultiCompartmentLocalBasisNode<Tree> ||
489 Concept::MultiMembraneLocalBasisNode<Tree>
490 void initialize_nodes(
const Tree& tree,
auto fname)
492 for (std::size_t c = 0; c != tree.degree(); ++c)
493 initialize_nodes(tree.child(c), fname, c);
496 void initialize_nodes(
const Concept::MultiCompartmentMembraneLocalBasisNode
auto& tree,
499 using namespace Indices;
500 initialize_nodes(tree.child(_0), fname);
501 initialize_nodes(tree.child(_1), fname);
504 void initialize_nodes(
const Concept::CompartmentMembraneLocalBasisNode
auto& tree,
auto fname)
506 using namespace Indices;
507 initialize_nodes(tree.child(_0), fname, 0);
508 initialize_nodes(tree.child(_1), fname, 0);
511 using ScalarTag = index_constant<0>;
512 using VectorTag = index_constant<1>;
513 using TensorApplyTag = index_constant<2>;
516 static fu2::unique_function<Scalar() const noexcept> make_functor(
517 const LocalEquations& local_eqs,
518 const FunctorFactory<dim>& functor_factory,
519 std::string_view prefix,
520 const ParameterTree& config,
524 return functor_factory.make_scalar(prefix, config, local_eqs, codim);
528 static fu2::unique_function<Vector() const noexcept> make_functor(
529 const LocalEquations& local_eqs,
530 const FunctorFactory<dim>& functor_factory,
531 std::string_view prefix,
532 const ParameterTree& config,
536 return functor_factory.make_vector(prefix, config, local_eqs, codim);
540 static fu2::unique_function<Vector(Vector)
const noexcept> make_functor(
541 const LocalEquations& local_eqs,
542 const FunctorFactory<dim>& functor_factory,
543 std::string_view prefix,
544 const ParameterTree& config,
548 return functor_factory.make_tensor_apply(prefix, config, local_eqs, codim);
552 template<
class Signature>
553 static auto set_differentiable_function(
const LocalEquations& local_eqs,
554 const FunctorFactory<dim>& functor_factory,
555 CompartmentDifferentiableFunction<Signature>& function,
556 std::string_view prefix,
557 const ParameterTree& function_config,
560 function = CompartmentDifferentiableFunction<Signature>{ make_functor(
561 local_eqs, functor_factory, prefix, function_config, 0, range_tag) };
562 std::set<std::string_view> debug_set;
563 const auto& jac_config = function_config.sub(
"jacobian");
564 for (
const auto& compartment_fncs_jac : local_eqs._nodes[Indices::_0])
565 for (
const auto& component_fncs_jac : compartment_fncs_jac)
566 if (jac_config.hasSub(component_fncs_jac.name))
568 make_functor(local_eqs,
570 fmt::format(
"{}.jacobian.{}", prefix, component_fncs_jac.name),
571 jac_config.sub(component_fncs_jac.name),
574 function.compartment_jacobian.emplace_back(std::move(jac), component_fncs_jac);
575 debug_set.insert(component_fncs_jac.name);
577 if (jac_config.getSubKeys().size() != debug_set.size())
578 spdlog::warn(
"Some sub-sections in \"jacobian\" section are being ignored");
582 template<
class Signature>
583 static auto set_differentiable_function(
const LocalEquations& local_eqs,
584 const FunctorFactory<dim>& functor_factory,
585 MembraneDifferentiableFunction<Signature>& function,
586 std::string_view prefix,
587 const ParameterTree& function_config,
590 function = MembraneDifferentiableFunction<Signature>{ make_functor(
591 local_eqs, functor_factory, prefix, function_config, 1, range_tag) };
592 const auto& jac_config = function_config.sub(
"jacobian");
593 for (
const auto& compartment_fncs_jac :local_eqs. _nodes[Indices::_0])
594 for (
const auto& component_fncs_jac : compartment_fncs_jac)
595 if (jac_config.hasSub(component_fncs_jac.name))
597 make_functor(local_eqs,
599 fmt::format(
"{}.jacobian.{}", prefix, component_fncs_jac.name),
600 jac_config.sub(component_fncs_jac.name),
603 function.compartment_jacobian.emplace_back(std::move(jac), component_fncs_jac);
604 for (
const auto& membrane_fncs_jac : local_eqs._nodes[Indices::_1])
605 for (
const auto& component_fncs_jac : membrane_fncs_jac)
606 if (jac_config.hasSub(component_fncs_jac.name))
608 make_functor(local_eqs,
610 fmt::format(
"{}.jacobian.{}", prefix, component_fncs_jac.name),
611 jac_config.sub(component_fncs_jac.name),
614 function.membrane_jacobian.emplace_back(std::move(jac), component_fncs_jac);
617 void configure_eqs(
const ParameterTree& eqs_cfg,
618 std::shared_ptr<
const FunctorFactory<dim>> functor_factory,
619 BitFlags<FactoryFalgs> opts)
622 _nodes, [&]<
class Node>(std::vector<std::vector<Node>>& compartments_fncs,
auto l) {
623 for (
auto& compartment_fncs : compartments_fncs) {
624 for (Node& component_fncs : compartment_fncs) {
626 const auto& component_config = eqs_cfg.sub(component_fncs.name,
true);
628 if (opts.test(FactoryFalgs::Diffusion) and component_config.hasSub(
"cross_diffusion")) {
629 const auto& cross_diffusion_config = component_config.sub(
"cross_diffusion");
630 std::set<std::string_view> debug_set;
631 for (const auto& compartment_fncs_cross_diff : compartments_fncs)
632 for (const auto& component_fncs_cross_diff : compartment_fncs_cross_diff)
633 if (cross_diffusion_config.hasSub(component_fncs_cross_diff.name)) {
634 const auto& diffusion_config =
635 cross_diffusion_config.sub(component_fncs_cross_diff.name);
636 debug_set.insert(component_fncs_cross_diff.name);
637 auto& cross_diffusion = component_fncs.cross_diffusion.emplace_back(
638 nullptr, component_fncs_cross_diff);
639 set_differentiable_function(
640 std::as_const(*this),
643 fmt::format(
"{}.cross_diffusion", component_fncs.name),
646 if (not cross_diffusion)
647 component_fncs.cross_diffusion.pop_back();
649 if (cross_diffusion_config.getSubKeys().size() != debug_set.size())
651 "Some sub-sections in \"{}.cross_diffusion\" section are being ignored",
652 component_fncs.name);
655 if (opts.test(FactoryFalgs::Reaction) and component_config.hasSub(
"reaction"))
656 set_differentiable_function(std::as_const(*
this),
658 component_fncs.reaction,
659 fmt::format(
"{}.reaction", component_fncs.name),
660 component_config.sub(
"reaction"),
663 if (opts.test(FactoryFalgs::Velocity) and component_config.hasSub(
"velocity"))
664 set_differentiable_function(std::as_const(*
this),
666 component_fncs.velocity,
667 fmt::format(
"{}.velocity", component_fncs.name),
668 component_config.sub(
"velocity"),
671 if (opts.test(FactoryFalgs::Storage) and component_config.hasSub(
"storage"))
672 set_differentiable_function(std::as_const(*
this),
674 component_fncs.storage,
675 fmt::format(
"{}.storage", component_fncs.name),
676 component_config.sub(
"storage"),
679 if (opts.test(FactoryFalgs::Outflow) and component_config.hasSub(
"outflow")) {
681 throw format_exception(NotImplemented{},
682 "Outflow for membranes is not implemented");
683 const auto& boundaries_config = component_config.sub(
"outflow");
684 for (std::size_t i = 0; i != _compartment_names[l].size(); ++i) {
685 if (boundaries_config.hasSub(_compartment_names[l][i])) {
686 component_fncs.outflow.resize(_compartment_names[l].size());
687 set_differentiable_function(
688 std::as_const(*
this),
690 component_fncs.outflow[i],
691 fmt::format(
"{}.outflow.{}", component_fncs.name, _compartment_names[l][i]),
692 boundaries_config.sub(_compartment_names[l][i]),
Bit flags for enumerators.
Definition: bit_flags.hh:47
static constexpr BitFlags< Enum > no_flags()
Bitflag with all flags turned off.
Definition: bit_flags.hh:70
constexpr bool any() const
Checks if any flags are set to true.
Definition: bit_flags.hh:207
Container for cell data of a grid view.
Definition: cell_data.hh:25
Definition: local_equations.hh:39
static std::unique_ptr< LocalEquations > make_stiffness(const PDELab::Concept::LocalBasis auto &lbasis, const ParameterTree &eqs_cfg, const ParameterTree &domain_cfg, std::shared_ptr< const FunctorFactory< dim > > functor_factory, std::shared_ptr< const CellData< CellDataGridView, CellDataType > > grid_cell_data)
Definition: local_equations.hh:323
void forEachValue(Codim< 1 >, const std::function< void(std::string_view, const FieldVector< double, 1 > &, const FieldVector< double, dim - 1 > &)> &) const override
Definition: local_equations.hh:393
LocalEquations & operator=(LocalEquations &&)=delete
void clear()
Definition: local_equations.hh:401
auto & get_gradient(const Tree &tree)
Definition: local_equations.hh:376
const auto & nodes() const
Definition: local_equations.hh:416
LocalEquations(const LocalEquations &)=delete
LocalEquations(LocalEquations &&)=delete
LocalEquations & operator=(const LocalEquations &)=delete
static std::unique_ptr< LocalEquations > make_mass(const PDELab::Concept::LocalBasis auto &lbasis, const ParameterTree &eqs_cfg, const ParameterTree &domain_cfg, std::shared_ptr< const FunctorFactory< dim > > functor_factory, std::shared_ptr< const CellData< CellDataGridView, CellDataType > > grid_cell_data)
Definition: local_equations.hh:341
CompartmentVectorFunction domain_deformation
Definition: local_equations.hh:418
void forEachValue(Codim< 0 >, const std::function< void(std::string_view, const FieldVector< double, 1 > &, const FieldVector< double, dim > &)> &apply) const override
Definition: local_equations.hh:381
virtual ~LocalEquations() override=default
static std::unique_ptr< LocalEquations > make(const LocalBasis &lbasis, const ParameterTree &eqs_cfg={}, const ParameterTree &domain_cfg={}, std::shared_ptr< const FunctorFactory< dim > > functor_factory=nullptr, std::shared_ptr< const CellData< CellDataGridView, CellDataType > > grid_cell_data=nullptr, BitFlags< FactoryFalgs > opts=BitFlags< FactoryFalgs >::no_flags())
Definition: local_equations.hh:264
void debug() const
Definition: local_equations.hh:407
auto & get_value(const Tree &tree)
Definition: local_equations.hh:368
const auto & get_equation(const Tree &tree) const
Definition: local_equations.hh:360
Definition: functor_factory.hh:24
Definition: functor_factory_parser.impl.hh:22
Definition: compartment_tree.hh:30
Definition: compartment_tree.hh:86
Definition: compartment_tree.hh:20
Definition: compartment_tree.hh:36
Definition: compartment_tree.hh:92
Definition: compartment_tree.hh:25
Definition: compartment_tree.hh:42
Definition: compartment_tree.hh:61
Definition: factory.hh:28
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition: exceptions.hh:23
Definition: local_equations.hh:127
std::vector< CompartmentPartialDerivative< Signature > > compartment_jacobian
Definition: local_equations.hh:128
Definition: local_equations.hh:150
const CompartmentNode & wrt
Definition: local_equations.hh:156
CompartmentDiffusionApply(fu2::unique_function< Vector(Vector) const noexcept > &&callable, const CompartmentNode &_wrt)
Definition: local_equations.hh:151
Definition: local_equations.hh:173
CompartmentScalarFunction storage
Definition: local_equations.hh:181
std::vector< MembraneScalarFunction > outflow
Definition: local_equations.hh:185
std::vector< CompartmentDiffusionApply > cross_diffusion
Definition: local_equations.hh:184
std::string name
Definition: local_equations.hh:178
void debug() const
Definition: local_equations.hh:195
CompartmentScalarFunction reaction
Definition: local_equations.hh:180
CompartmentNode(Scalar &value, Vector &gradient, CompartmentPath path, const std::string &name)
Definition: local_equations.hh:187
Vector & gradient
Definition: local_equations.hh:175
const Concept::CompartmentScalarLocalBasisNode auto & to_local_basis_node(const PDELab::Concept::LocalBasis auto &lbasis) const
Definition: local_equations.hh:189
CompartmentPath path
Definition: local_equations.hh:177
Scalar & value
Definition: local_equations.hh:174
CompartmentVectorFunction velocity
Definition: local_equations.hh:182
Definition: local_equations.hh:103
const CompartmentNode & wrt
Definition: local_equations.hh:110
CompartmentPartialDerivative(fu2::unique_function< Signature > &&callable, const CompartmentNode &_wrt)
Definition: local_equations.hh:104
Definition: local_equations.hh:133
std::vector< MembranePartialDerivative< Signature > > membrane_jacobian
Definition: local_equations.hh:135
std::vector< CompartmentPartialDerivative< Signature > > compartment_jacobian
Definition: local_equations.hh:134
Definition: local_equations.hh:161
const MembraneNode & wrt
Definition: local_equations.hh:169
MembraneDiffusionApply(fu2::unique_function< Vector(Vector) const noexcept > &&callable, const MembraneNode &_wrt)
Definition: local_equations.hh:162
Definition: local_equations.hh:219
void debug() const
Definition: local_equations.hh:240
MembranePath path
Definition: local_equations.hh:223
Scalar & value
Definition: local_equations.hh:220
std::vector< MembraneDiffusionApply > cross_diffusion
Definition: local_equations.hh:231
bool is_linear
Definition: local_equations.hh:225
MembraneScalarFunction storage
Definition: local_equations.hh:228
std::string name
Definition: local_equations.hh:224
Vector & gradient
Definition: local_equations.hh:221
const Concept::MembraneSubEntitiesLocalBasisNode auto & to_local_basis_node(const PDELab::Concept::LocalBasis auto &lbasis) const
Definition: local_equations.hh:234
MembraneScalarFunction reaction
Definition: local_equations.hh:227
std::vector< MembraneScalarFunction > outflow
Definition: local_equations.hh:232
MembraneVectorFunction velocity
Definition: local_equations.hh:229
Definition: local_equations.hh:115
const MembraneNode & wrt
Definition: local_equations.hh:122
MembranePartialDerivative(fu2::unique_function< Signature > &&callable, const MembraneNode &_wrt)
Definition: local_equations.hh:116
Definition: local_domain.hh:15