Dune::Copasi
Loading...
Searching...
No Matches
make_initial.hh
Go to the documentation of this file.
1#ifndef DUNE_COPASI_MODEL_MAKE_INITIAL_HH
2#define DUNE_COPASI_MODEL_MAKE_INITIAL_HH
3
7
8#include <dune/pdelab/common/trace.hh>
9
10#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
11
12#include <dune/grid/concepts/grid.hh>
13
14#include <dune/common/exceptions.hh>
15
16#include <spdlog/spdlog.h>
17
18#include <optional>
19#include <string>
20#include <unordered_map>
21
22namespace Dune::Copasi {
23
24template<class GridFunction, Dune::Concept::Grid Grid>
25[[nodiscard]] inline static std::unordered_map<std::string, GridFunction>
26make_initial(const Grid& grid,
27 const ParameterTree& config,
29{
30 TRACE_EVENT("dune", "InitialCondition");
31
32 std::unordered_map<std::string, GridFunction> functions;
33
34 auto time = config.get("time_step_operator.time_begin", double{ 0. });
35
36 const auto& compartments_config = config.sub("compartments", true);
37
38 for (const auto& compartment : compartments_config.getSubKeys()) {
39
40 // find grid view to use
41 auto sub_grid_view = [&]() {
43 if constexpr (Concept::MultiDomainGrid<Grid>) {
44 using SubDomainIndex = typename Grid::SubDomainIndex;
45 auto const domain_id = compartment_config.template get<SubDomainIndex>("id");
46 return grid.subDomain(domain_id).leafGridView();
47 } else {
48 auto const domain_id = compartment_config.template get<std::size_t>("id");
49 if (domain_id != 0) {
50 throw format_exception(IOError{}, "Compartment ID does not exist in grid");
51 }
52 return grid.leafGridView();
53 }
54 }();
55
56 auto components = config.sub("scalar_field", true).getSubKeys();
57 for (const auto& component : components) {
58 auto prefix = fmt::format("scalar_field.{}.initial", component);
59 if (not config.hasSub(prefix)) {
60 continue;
61 }
62 auto local_domain = std::make_shared<LocalDomain<Grid::dimensionworld>>();
63 local_domain->time = time;
64 auto fcn = functor_factory.make_scalar(prefix, config.sub(prefix), *local_domain);
65 if (not fcn) {
66 continue;
67 }
68 auto shared_fnc = std::make_shared<decltype(fcn)>(std::move(fcn));
69 using GridDomain =
70 typename decltype(sub_grid_view)::template Codim<0>::Geometry::GlobalCoordinate;
71 std::function<double(GridDomain)> analytic_fcn =
72 [lcl_domain = std::move(local_domain),
74 lcl_domain->position = global_position;
75 return std::invoke(*_fnc);
76 };
77
78 auto grid_function =
79 Dune::Functions::makeAnalyticGridViewFunction(std::move(analytic_fcn), sub_grid_view);
80
81 std::string const component_id = fmt::format("{}.{}", compartment, component);
82 auto [_, inserted] = functions.try_emplace(component_id, std::move(grid_function));
83 if (not inserted) {
84 throw format_exception(IOError{}, "Two or more components share the same name");
85 }
86 }
87 }
88
89 return functions;
90}
91
92} // namespace Dune::Copasi
93
94#endif // DUNE_COPASI_MODEL_MAKE_INITIAL_HH
Definition axis_names.hh:7
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition exceptions.hh:23
constexpr bool is_bitflags_v
Alias for Bitflag indicator.
Definition bit_flags.hh:24