Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
symengine_scalar_operations.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#include <deal.II/base/config.h>
17
18#ifdef DEAL_II_WITH_SYMENGINE
19
24
25# include <symengine/real_double.h>
26
28
29namespace Differentiation
30{
31 namespace SD
32 {
33 namespace SE = ::SymEngine;
34
35
36 /* ------------------------- Symbol creation -----------------------*/
37
38
39 Expression
40 make_symbol(const std::string &symbol)
41 {
42 return Expression(symbol);
43 }
44
45
46# ifndef DOXYGEN
47 Expression
48 make_symbolic_function(const std::string &symbol,
49 const SD::types::symbol_vector &arguments)
50 {
51 return Expression(symbol, arguments);
52 }
53
54
55 Expression
56 make_symbolic_function(const std::string &symbol,
57 const SD::types::substitution_map &arguments)
58 {
59 return make_symbolic_function(symbol,
61 }
62# endif
63
64
65 /* --------------------------- Differentiation ------------------------- */
66
67
68 Expression
69 differentiate(const Expression &func, const Expression &op)
70 {
71 return func.differentiate(op);
72 }
73
74
75 /* ---------------- Symbol map creation and manipulation --------------*/
76
77
78# ifndef DOXYGEN
79 namespace internal
80 {
81 bool
82 is_valid_substitution_symbol(const SE::Basic &entry)
83 {
84 // Allow substitution of a symbol
85 // It is pretty clear as to why this is wanted...
86 if (SE::is_a<SE::Symbol>(entry))
87 return true;
88
89 // Allow substitution of a function symbol
90 // If desired, we can transform general but undefined functional
91 // relationships to an explicit form that is concrete. This is
92 // required for a symbolic expression to be parsed by a Lambda or LLVM
93 // optimizer.
94 if (SE::is_a<SE::FunctionSymbol>(entry))
95 return true;
96
97 // Allow substitution of the explicit expression of the derivative of
98 // an implicitly defined symbol (i.e. the result of the derivative of
99 // a FunctionSymbol).
100 if (SE::is_a<SE::Derivative>(entry))
101 return true;
102
103 // When performing tensor differentiation, one may end up with a
104 // coefficient of one half due to symmetry operations, e.g.
105 // 0.5*Derivative(Qi_00(C_11, C_00, C_01), C_01)
106 // So we explicitly check for this exact case
107 if (SE::is_a<SE::Mul>(entry))
108 {
109 const SE::Mul &entry_mul = SE::down_cast<const SE::Mul &>(entry);
110 // Check that the factor is a half...
111 if (SE::eq(*(entry_mul.get_coef()), *SE::real_double(0.5)))
112 {
113 // ...and that there is only one entry and that its a
114 // Derivative type
115 const SE::map_basic_basic &entry_mul_dict =
116 entry_mul.get_dict();
117 if (entry_mul_dict.size() == 1 &&
118 SE::is_a<SE::Derivative>(*(entry_mul_dict.begin()->first)))
119 return true;
120 }
121 }
122
123 return false;
124 }
125
126
127 void
129 types::substitution_map &substitution_map,
130 const SymEngine::RCP<const SymEngine::Basic> &symbol,
131 const SymEngine::RCP<const SymEngine::Basic> &value)
132 {
133 Assert(
136 "Substitution with a number that does not represent a symbol or symbolic derivative"));
137
138 auto it_sym = substitution_map.find(Expression(symbol));
139 Assert(it_sym != substitution_map.end(),
140 ExcMessage("Did not find this symbol in the map."));
141
142 it_sym->second = Expression(value);
143 }
144
145 } // namespace internal
146
147
148 void
150 const Expression &symbol,
151 const Expression &value)
152 {
153 internal::set_value_in_symbol_map(substitution_map,
154 symbol.get_RCP(),
155 value.get_RCP());
156 }
157
158
159 void
161 const types::substitution_map &symbol_values)
162 {
163 for (const auto &entry : symbol_values)
165 }
166# endif
167
168
169 /* ------------------ Symbol substitution map creation ----------------*/
170
171
173 make_substitution_map(const Expression &symbol, const Expression &value)
174 {
175 types::substitution_map substitution_map;
176 add_to_substitution_map(substitution_map, symbol, value);
177 return substitution_map;
178 }
179
180
181# ifndef DOXYGEN
182 /* ---------------- Symbolic substitution map enlargement --------------*/
183
184
185 void
187 const types::substitution_map &symb_map_in)
188 {
189 // Do this by hand so that we can perform some sanity checks
190 for (const auto &entry : symb_map_in)
191 {
192 const typename types::substitution_map::const_iterator it_other =
193 symb_map_out.find(entry.first);
194 if (it_other == symb_map_out.end())
195 symb_map_out.insert(std::make_pair(entry.first, entry.second));
196 else
197 {
198 Assert(SE::eq(*(entry.second.get_RCP()),
199 *(it_other->second.get_RCP())),
200 ExcMessage("Key already in map, but values don't match"));
201 }
202 }
203 }
204
205
206 /* ---------------- Symbol substitution and evaluation --------------*/
207
208
211 const bool force_cyclic_dependency_resolution)
212 {
213 types::substitution_map symbol_values_resolved = symbol_values;
214 const std::size_t size = symbol_values.size();
215 (void)size;
216 for (auto &entry : symbol_values_resolved)
217 {
218 // Perform dictionary-based substitution to
219 // resolve all explicit relations in a map.
220 // Instead of checking by value (and thus having
221 // to store a temporary value), we check to see
222 // if the hash of the map entry changes.
223 Expression &out = entry.second;
224 SE::hash_t hash_old;
225 SE::hash_t hash_new = out.get_RCP()->hash();
226 unsigned int iter = 0;
227 do
228 {
229 // Write the substituted value straight back
230 // into the map.
231 if (force_cyclic_dependency_resolution)
232 {
233 // Here we substitute in symbol_values_resolved instead of
234 // symbol_values, in order to break any cyclic dependencies.
235 // The earlier entries in the dictionary are in this way
236 // guaranteed to be resolved before any subsequent entries,
237 // thereby breaking the dependency loop.
238 out = out.substitute(symbol_values_resolved);
239 }
240 else
241 {
242 out = out.substitute(symbol_values);
243 }
244
245 // Compute and store the hash of the new object
246 hash_old = hash_new;
247 hash_new = out.get_RCP()->hash();
249 iter < size,
251 "Unresolvable cyclic dependency detected in substitution map."));
252 ++iter;
253 }
254 while (hash_new != hash_old);
255 }
256
257 return symbol_values_resolved;
258 }
259
260
261 Expression
262 substitute(const Expression &expression,
263 const types::substitution_map &substitution_map)
264 {
265 return expression.substitute(substitution_map);
266 }
267# endif // DOXYGEN
268
269 } // namespace SD
270} // namespace Differentiation
271
273
274#endif // DEAL_II_WITH_SYMENGINE
Expression substitute(const types::substitution_map &substitution_values) const
Expression differentiate(const Expression &symbol) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4614
Point< 2 > first
Definition grid_out.cc:4613
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
SD::types::symbol_vector extract_symbols(const SD::types::substitution_map &substitution_values)
bool is_valid_substitution_symbol(const SymEngine::Basic &entry)
void set_value_in_symbol_map(types::substitution_map &substitution_map, const SymEngine::RCP< const SymEngine::Basic > &symbol, const SymEngine::RCP< const SymEngine::Basic > &value)
std::vector< SD::Expression > symbol_vector
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
void merge_substitution_maps(types::substitution_map &substitution_map_out, const types::substitution_map &substitution_map_in)
Expression differentiate(const Expression &f, const Expression &x)
Expression make_symbolic_function(const std::string &symbol, const types::symbol_vector &arguments)
void set_value_in_symbol_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
void add_to_substitution_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
types::substitution_map resolve_explicit_dependencies(const types::substitution_map &substitution_map, const bool force_cyclic_dependency_resolution=false)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
Expression make_symbol(const std::string &symbol)
const InputIterator OutputIterator out
Definition parallel.h:167