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