Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
symengine_scalar_operations.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 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 
21 # include <deal.II/differentiation/sd/symengine_number_types.h>
22 # include <deal.II/differentiation/sd/symengine_scalar_operations.h>
23 # include <deal.II/differentiation/sd/symengine_types.h>
24 # include <deal.II/differentiation/sd/symengine_utilities.h>
25 
26 # include <symengine/real_double.h>
27 
28 DEAL_II_NAMESPACE_OPEN
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  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,
60  SD::Utilities::extract_symbols(arguments));
61  }
62 
63 
64  /* --------------------------- Differentiation ------------------------- */
65 
66 
67  Expression
68  differentiate(const Expression &func, const Expression &op)
69  {
70  return func.differentiate(op);
71  }
72 
73 
74  /* ---------------- Symbol map creation and manipulation --------------*/
75 
76 
77  namespace internal
78  {
79  bool
80  is_valid_substitution_symbol(const SE::Basic &entry)
81  {
82  // Allow substitution of a symbol
83  // It is pretty clear as to why this is wanted...
84  if (SE::is_a<SE::Symbol>(entry))
85  return true;
86 
87  // Allow substitution of a function symbol
88  // If desired, we can transform general but undefined functional
89  // relationships to an explicit form that is concrete. This is
90  // required for a symbolic expression to be parsed by a Lambda or LLVM
91  // optimizer.
92  if (SE::is_a<SE::FunctionSymbol>(entry))
93  return true;
94 
95  // Allow substitution of the explicit expression of the derivative of
96  // an implicitly defined symbol (i.e. the result of the derivative of
97  // a FunctionSymbol).
98  if (SE::is_a<SE::Derivative>(entry))
99  return true;
100 
101  // When performing tensor differentiation, one may end up with a
102  // coefficient of one half due to symmetry operations, e.g.
103  // 0.5*Derivative(Qi_00(C_11, C_00, C_01), C_01)
104  // So we explicitly check for this exact case
105  if (SE::is_a<SE::Mul>(entry))
106  {
107  const SE::Mul &entry_mul = SE::down_cast<const SE::Mul &>(entry);
108  // Check that the factor is a half...
109  if (SE::eq(*(entry_mul.get_coef()), *SE::real_double(0.5)))
110  {
111  // ...and that there is only one entry and that its a
112  // Derivative type
113  const SE::map_basic_basic &entry_mul_dict =
114  entry_mul.get_dict();
115  if (entry_mul_dict.size() == 1 &&
116  SE::is_a<SE::Derivative>(*(entry_mul_dict.begin()->first)))
117  return true;
118  }
119  }
120 
121  return false;
122  }
123 
124 
125  void
127  types::substitution_map & substitution_map,
128  const SymEngine::RCP<const SymEngine::Basic> &symbol,
129  const SymEngine::RCP<const SymEngine::Basic> &value)
130  {
131  Assert(
132  internal::is_valid_substitution_symbol(*symbol),
133  ExcMessage(
134  "Substitution with a number that does not represent a symbol or symbolic derivative"));
135 
136  auto it_sym = substitution_map.find(Expression(symbol));
137  Assert(it_sym != substitution_map.end(),
138  ExcMessage("Did not find this symbol in the map."));
139 
140  it_sym->second = Expression(value);
141  }
142 
143  } // namespace internal
144 
145 
146  void
147  set_value_in_symbol_map(types::substitution_map &substitution_map,
148  const Expression & symbol,
149  const Expression & value)
150  {
151  internal::set_value_in_symbol_map(substitution_map,
152  symbol.get_RCP(),
153  value.get_RCP());
154  }
155 
156 
157  void
158  set_value_in_symbol_map(types::substitution_map & substitution_map,
159  const types::substitution_map &symbol_values)
160  {
161  for (const auto &entry : symbol_values)
162  set_value_in_symbol_map(substitution_map, entry.first, entry.second);
163  }
164 
165 
166  /* ------------------ Symbol substitution map creation ----------------*/
167 
168 
169  types::substitution_map
170  make_substitution_map(const Expression &symbol, const Expression &value)
171  {
172  types::substitution_map substitution_map;
173  add_to_substitution_map(substitution_map, symbol, value);
174  return substitution_map;
175  }
176 
177 
178  /* ---------------- Symbolic substitution map enlargement --------------*/
179 
180 
181  void
182  merge_substitution_maps(types::substitution_map & symb_map_out,
183  const types::substitution_map &symb_map_in)
184  {
185  // Do this by hand so that we can perform some sanity checks
186  for (const auto &entry : symb_map_in)
187  {
188  const typename types::substitution_map::const_iterator it_other =
189  symb_map_out.find(entry.first);
190  if (it_other == symb_map_out.end())
191  symb_map_out.insert(std::make_pair(entry.first, entry.second));
192  else
193  {
194  Assert(SE::eq(*(entry.second.get_RCP()),
195  *(it_other->second.get_RCP())),
196  ExcMessage("Key already in map, but values don't match"));
197  }
198  }
199  }
200 
201 
202  /* ---------------- Symbol substitution and evaluation --------------*/
203 
204 
205  types::substitution_map
206  resolve_explicit_dependencies(const types::substitution_map &symbol_values,
207  const bool force_cyclic_dependency_resolution)
208  {
209  types::substitution_map symbol_values_resolved = symbol_values;
210  const std::size_t size = symbol_values.size();
211  (void)size;
212  for (auto &entry : symbol_values_resolved)
213  {
214  // Perform dictionary-based substitution to
215  // resolve all explicit relations in a map.
216  // Instead of checking by value (and thus having
217  // to store a temporary value), we check to see
218  // if the hash of the map entry changes.
219  Expression & out = entry.second;
220  SE::hash_t hash_old;
221  SE::hash_t hash_new = out.get_RCP()->hash();
222  unsigned int iter = 0;
223  do
224  {
225  // Write the substituted value straight back
226  // into the map.
227  if (force_cyclic_dependency_resolution)
228  {
229  // Here we substitute in symbol_values_resolved instead of
230  // symbol_values, in order to break any cyclic dependencies.
231  // The earlier entries in the dictionary are in this way
232  // guaranteed to be resolved before any subsequent entries,
233  // thereby breaking the dependency loop.
234  out = out.substitute(symbol_values_resolved);
235  }
236  else
237  {
238  out = out.substitute(symbol_values);
239  }
240 
241  // Compute and store the hash of the new object
242  hash_old = std::move(hash_new);
243  hash_new = out.get_RCP()->hash();
244  AssertThrow(
245  iter < size,
246  ExcMessage(
247  "Unresolvable cyclic dependency detected in substitution map."));
248  ++iter;
249  }
250  while (hash_new != hash_old);
251  }
252 
253  return symbol_values_resolved;
254  }
255 
256 
257  Expression
258  substitute(const Expression & expression,
259  const types::substitution_map &substitution_map)
260  {
261  return expression.substitute(substitution_map);
262  }
263 
264  } // namespace SD
265 } // namespace Differentiation
266 
267 DEAL_II_NAMESPACE_CLOSE
268 
269 #endif // DEAL_II_WITH_SYMENGINE
Expression differentiate(const Expression &symbol) const
void set_value_in_symbol_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
void merge_substitution_maps(types::substitution_map &substitution_map_out, const types::substitution_map &substitution_map_in)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
Expression make_symbolic_function(const std::string &symbol, const types::symbol_vector &arguments)
types::substitution_map resolve_explicit_dependencies(const types::substitution_map &substitution_map, const bool force_cyclic_dependency_resolution=false)
Expression differentiate(const Expression &f, const Expression &x)
Expression substitute(const types::substitution_map &substitution_values) const
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
void add_to_substitution_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
Expression make_symbol(const std::string &symbol)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)