1#ifndef DUNE_COPASI_MODEL_DIFFUSION_REACTION_SINGLE_COMPARTMENT_IMPL_HH
2#define DUNE_COPASI_MODEL_DIFFUSION_REACTION_SINGLE_COMPARTMENT_IMPL_HH
17#include <dune/pdelab/basis/backend/istl.hh>
18#include <dune/pdelab/basis/basis.hh>
19#include <dune/pdelab/basis/discrete_global_function.hh>
20#include <dune/pdelab/common/partition/identity.hh>
21#include <dune/pdelab/common/trace.hh>
22#include <dune/pdelab/concepts/multiindex.hh>
24#include <dune/grid/io/file/vtk.hh>
26#include <spdlog/spdlog.h>
28#ifdef DUNE_COPASI_PRECOMPILED_MODE
29#warning "Including this file in pre-compiled mode may defeat the purpose of pre-compilation"
36ModelSingleCompartment<Traits>::get_entity_set(
const Grid& grid, std::size_t subdomain)
37 -> CompartmentEntitySet
39 if constexpr (std::same_as<typename Grid::LeafGridView, CompartmentEntitySet>) {
40 return grid.leafGridView();
41 }
else if constexpr (Concept::SubDomainGrid<typename CompartmentEntitySet::Grid>) {
42 static_assert(std::same_as<typename Grid::SubDomainGrid::LeafGridView, CompartmentEntitySet>);
43 return grid.subDomain(subdomain).leafGridView();
45 throw format_exception(NotImplemented{},
"Not known mapping from Grid to CompartmentEntitySet");
52 const std::unordered_map<std::string, GridFunction>& initial)
const
54 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>,
CompartmentPreBasis>;
55 using CoefficientsBackend = PDELab::ISTLUniformBackend<ScalarQuantity>;
56 using Coefficients =
typename CompartmentBasis::template Container<CoefficientsBackend>;
57 const auto& basis = any_cast<const CompartmentBasis&>(state.basis);
58 auto& coefficients = any_cast<Coefficients&>(state.coefficients);
60 Dune::Copasi::interpolate(basis, coefficients, initial);
66 const CompartmentEntitySet& entity_set,
67 std::string_view name,
68 const ParameterTree& scalar_field_config,
71 spdlog::info(
"Setup basis functions for component '{}'", name);
72 auto scalar_field_pre_basis =
73 ScalarPreBasis{ ScalarMergingStrategy{ entity_set },
74 std::make_shared<ScalarFiniteElementMap>(entity_set),
75 Constraints{boundary_mapper, scalar_field_config.sub(
"constrain"), functor_factory} };
76 scalar_field_pre_basis.name(name);
77 return scalar_field_pre_basis;
84 std::string_view compartment_name,
85 const std::vector<std::string>& scalar_field_names,
86 const ParameterTree& scalar_fields_config,
89 spdlog::info(
"Setup compartment basis functions for compartment '{}'", compartment_name);
91 if (entity_set.size(0) == 0)
92 spdlog::warn(
"Compartment '{}' is empty", compartment_name);
94 auto boundary_mapper = std::make_shared<BoundaryEntityMapper<CompartmentEntitySet>>(entity_set);
95 std::vector<ScalarPreBasis> scalar_field_pre_basis;
96 for (
const auto& name : scalar_field_names) {
97 scalar_field_pre_basis.push_back(make_scalar_field_pre_basis(boundary_mapper, entity_set, name, scalar_fields_config.sub(name), functor_factory));
102 compartment_pre_basis.name(compartment_name);
104 spdlog::info(
"No. of components on '{}' compartment: {}",
105 compartment_pre_basis.name(),
106 compartment_pre_basis.degree());
108 return compartment_pre_basis;
111template<
class Traits>
115 const ParameterTree& config,
118 TRACE_EVENT(
"dune",
"Basis::SetUp");
119 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>, CompartmentPreBasis>;
120 const auto& compartments_config = config.sub(
"compartments",
true);
121 const auto& compartments = compartments_config.getSubKeys();
122 if (compartments.size() != 1) {
123 throw format_exception(InvalidStateException{},
"Config file should only have one compartment");
125 auto compartment = compartments.front();
126 auto entity_set = get_entity_set(
127 grid, compartments_config.sub(compartment,
true).template get<std::size_t>(
"id"));
129 std::vector<std::string> components;
130 const auto& scalar_fields_config = config.sub(
"scalar_field",
true);
131 for (
const auto& scalar_field : scalar_fields_config.getSubKeys()) {
132 if (scalar_fields_config.sub(scalar_field)[
"compartment"] == compartment) {
133 components.push_back(scalar_field);
136 auto comp_space = make_compartment_pre_basis(entity_set, compartment, components, scalar_fields_config, functor_factory);
137 state.basis = CompartmentBasis{ makeBasis(entity_set, comp_space) };
140template<
class Traits>
142ModelSingleCompartment<Traits>::setup_coefficient_vector(State& state)
144 spdlog::info(
"Setup coefficient vector");
145 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>, CompartmentPreBasis>;
146 using CoefficientsBackend = PDELab::ISTLUniformBackend<ScalarQuantity>;
147 using Coefficients =
typename CompartmentBasis::template Container<CoefficientsBackend>;
148 const auto& basis = any_cast<const CompartmentBasis&>(state.basis);
149 state.coefficients = Coefficients{ basis.makeContainer(CoefficientsBackend{}) };
152template<
class Traits>
155 const ParameterTree& config)
const
156 -> std::unique_ptr<State>
158 auto state_ptr = std::make_unique<State>();
159 setup_basis(*state_ptr, *grid, config, _functor_factory);
160 setup_coefficient_vector(*state_ptr);
161 state_ptr->grid = grid;
166template<
class Traits>
169 std::string_view name)
const
171 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>,
CompartmentPreBasis>;
172 using CoefficientsBackend = PDELab::ISTLUniformBackend<ScalarQuantity>;
173 using Coefficients =
typename CompartmentBasis::template Container<CoefficientsBackend>;
174 const auto& basis = any_cast<const CompartmentBasis&>(state->basis);
175 const auto& coeff = any_cast<const Coefficients&>(state->coefficients);
176 std::shared_ptr<const Coefficients>
const coeff_ptr(state, &coeff);
178 for (std::size_t component = 0; component != basis.degree(); ++component) {
179 auto leaf_space = basis.subSpace(TypeTree::treePath(std::size_t{ component }));
180 if (leaf_space.name() == name) {
181 return makeDiscreteGlobalBasisFunction(leaf_space, coeff_ptr);
184 throw format_exception(RangeError{},
"State doesn't contain any function with name: {}", name);
187template<
class Traits>
190 -> std::unordered_map<std::string, GridFunction>
192 return Dune::Copasi::make_initial<GridFunction>(grid, config, *_functor_factory);
195template<
class Traits>
198 const ParameterTree& config)
const
199 -> std::unique_ptr<PDELab::OneStep<State>>
201 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>,
CompartmentPreBasis>;
202 using CoefficientsBackend = PDELab::ISTLUniformBackend<ScalarQuantity>;
203 using Coefficients =
typename CompartmentBasis::template Container<CoefficientsBackend>;
204 using ResidualBackend = PDELab::ISTLUniformBackend<ResidualQuantity>;
205 using Residual =
typename CompartmentBasis::template Container<ResidualBackend>;
206 using LocalBasisTraits =
typename ScalarFiniteElementMap::Traits::FiniteElement::Traits::LocalBasisType::Traits;
208 const auto& basis = any_cast<const CompartmentBasis&>(state.basis);
210 std::shared_ptr one_step = DiffusionReaction::make_step_operator<LocalBasisTraits, Coefficients, Residual, ResidualQuantity, TimeQuantity>(config, basis, 1, _functor_factory, _cell_data);
213 auto type_erased_one_step = std::make_unique<PDELab::OperatorAdapter<State, State>>(
214 [one_step](PDELab::Operator<State, State>& self,
State& domain,
State& range)
mutable {
217 static_cast<PDELab::PropertyTree&
>(*one_step) =
218 static_cast<const PDELab::PropertyTree&
>(self);
219 auto& domain_coeff = any_cast<Coefficients&>(domain.coefficients);
220 auto& range_coeff = any_cast<Residual&>(range.coefficients);
221 return one_step->apply(domain_coeff, range_coeff);
225 static_cast<PDELab::PropertyTree&
>(*type_erased_one_step) =
226 static_cast<const PDELab::PropertyTree&
>(*one_step);
228 auto residual_ptr = std::make_shared<Residual>(basis.makeContainer(ResidualBackend{}));
229 type_erased_one_step->get(
"initial_residual") = residual_ptr;
230 type_erased_one_step->get(
"time") = state.time;
231 return type_erased_one_step;
234template<
class Traits>
237 -> std::map<std::string, double>
239 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>,
CompartmentPreBasis>;
240 using CoefficientsBackend = PDELab::ISTLUniformBackend<ScalarQuantity>;
241 using Coefficients =
typename CompartmentBasis::template Container<CoefficientsBackend>;
244 const auto& basis = any_cast<const CompartmentBasis&>(state.basis);
245 const auto& coeff = any_cast<const Coefficients&>(state.coefficients);
247 return Dune::Copasi::DiffusionReaction::reduce(PDELab::Execution::seq, basis, coeff, state.time, reduce_cfg, domain_cfg, _functor_factory);
250template<
class Traits>
253 const std::filesystem::path& path,
256 using CompartmentBasis = PDELab::Basis<PDELab::EntitySetPartitioner::Identity<CompartmentEntitySet>,
CompartmentPreBasis>;
257 using CoefficientsBackend = PDELab::ISTLUniformBackend<ScalarQuantity>;
258 using Coefficients =
typename CompartmentBasis::template Container<CoefficientsBackend>;
259 const auto& basis = any_cast<const CompartmentBasis&>(state.basis);
260 const auto& coeff = any_cast<const Coefficients&>(state.coefficients);
264 std::shared_ptr<const Coefficients>
const coeff_ptr = Dune::stackobject_to_shared_ptr(coeff);
267 auto path_entry = std::filesystem::directory_entry{ path };
268 if (not path_entry.exists()) {
269 spdlog::info(
"Creating output directory '{}'", path_entry.path().string());
270 std::error_code ec{ 0, std::generic_category() };
271 std::filesystem::create_directories(path_entry.path(), ec);
274 "\n Category: {}\nValue: {}\nMessage: {}\n",
275 ec.category().name(),
282 auto& timesteps = _writer_timesteps[path.string()];
283 std::string name = fmt::format(
"{}-{}", path.filename().string(), basis.name());
286 spdlog::info(
"Creating a time sequence file: '{}.pvd'", name);
288 spdlog::info(
"Overriding time sequence file: '{}.pvd'", name);
293 std::make_shared<VTKWriter<CompartmentEntitySet>>(basis.entitySet(), Dune::VTK::conforming);
294 auto sequential_writer =
295 VTKSequenceWriter<CompartmentEntitySet>{ writer, name, path.string(), path.string() };
296 sequential_writer.setTimeSteps(timesteps);
298 for (std::size_t component = 0; component != basis.degree(); ++component) {
299 auto const component_basis =
300 basis.subSpace(TypeTree::treePath(std::size_t{ component }));
301 sequential_writer.vtkWriter()->addVertexData(
302 makeDiscreteGlobalBasisFunction(component_basis, coeff_ptr),
303 VTK::FieldInfo{ component_basis.name(), VTK::FieldInfo::Type::scalar, 1 });
306 spdlog::info(
"Writing solution for {:.2f}s time stamp", state.time);
307 spdlog::info(
"Writing vtu file: '{0}/{0}-{1:0>5}.vtu'", name, timesteps.size());
308 TRACE_EVENT(
"dune",
"WriteVTK");
309 sequential_writer.write(state.time, Dune::VTK::base64);
310 sequential_writer.vtkWriter()->clear();
311 timesteps = sequential_writer.getTimeSteps();
313 if (coeff_ptr.use_count() != 1) {
315 InvalidStateException{},
316 "Fake shared pointer from coefficient vector may have been leaked outside of this"
Mapper of boundary entities.
Definition: boundary_entity_mapper.hh:26
Constraints parser and operator for a leaf basis.
Definition: constraints.hh:29
Definition: model_single_compartment.hh:26
std::unique_ptr< State > make_state(const std::shared_ptr< const Grid > &, const ParameterTree &) const override
Definition: model_single_compartment.impl.hh:154
typename Traits::CompartmentEntitySet CompartmentEntitySet
Definition: model_single_compartment.hh:37
typename Traits::CompartmentMergingStrategy CompartmentMergingStrategy
Definition: model_single_compartment.hh:42
static CompartmentPreBasis make_compartment_pre_basis(const CompartmentEntitySet &, std::string_view, const std::vector< std::string > &, const ParameterTree &={}, std::shared_ptr< const FunctorFactory< Grid::dimensionworld > >=nullptr)
Definition: model_single_compartment.impl.hh:82
typename Traits::Grid Grid
Definition: model_single_compartment.hh:34
typename Base::State State
Definition: model_single_compartment.hh:33
GridFunction make_compartment_function(const std::shared_ptr< const State > &, std::string_view) const override
Definition: model_single_compartment.impl.hh:168
std::map< std::string, double > reduce(const State &, const ParameterTree &, const ParameterTree &={}) const override
Definition: model_single_compartment.impl.hh:236
PDELab::PreBasisVector< CompartmentMergingStrategy, ScalarPreBasis > CompartmentPreBasis
Definition: model_single_compartment.hh:43
typename Base::GridFunction GridFunction
Definition: model_single_compartment.hh:46
void interpolate(State &, const std::unordered_map< std::string, GridFunction > &) const override
Definition: model_single_compartment.impl.hh:50
void write_vtk(const State &, const std::filesystem::path &, bool=true) const override
Definition: model_single_compartment.impl.hh:252
typename Traits::TimeQuantity TimeQuantity
Definition: model_single_compartment.hh:35
std::unordered_map< std::string, GridFunction > make_initial(const Grid &, const ParameterTree &) const override
Definition: model_single_compartment.impl.hh:189
std::unique_ptr< PDELab::OneStep< State > > make_step_operator(const State &, const ParameterTree &) const override
Definition: model_single_compartment.impl.hh:197
Definition: functor_factory.hh:24
Definition: factory.hh:28
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition: exceptions.hh:23
auto ostream2spdlog()
Standard output redirection to spdlog.
Definition: ostream_to_spdlog.hh:20