VertexCFD  0.0-dev
VertexCFD_Utils_SmoothMax.hpp
1 #ifndef VERTEXCFD_UTILS_SMOOTHMAX_HPP
2 #define VERTEXCFD_UTILS_SMOOTHMAX_HPP
3 
4 #include <Sacado.hpp>
5 #include <Sacado_Fad_Exp_Expression.hpp>
6 #include <Sacado_Fad_Exp_ExpressionTraits.hpp>
7 #include <Sacado_mpl_disable_if.hpp>
8 #include <Sacado_mpl_enable_if.hpp>
9 
10 namespace VertexCFD
11 {
12 namespace SmoothMath
13 {
14 //---------------------------------------------------------------------------//
15 // Implementation of smooth "max" using Sacado expressions.
16 // The third template parameter specifies whether the second parameter is a
17 // constant (non-FAD). The "max" function is commutative, so we do not need
18 // separate implementations of the expression class when T1 is FAD and T2 is
19 // constant AND when T1 is constant and T2 is FAD. If exactly one parameter is
20 // a FAD, it should be the first parameter. The reverse ordering (first
21 // parameter is constant) is handled within the specializations of the "max"
22 // function below.
23 //---------------------------------------------------------------------------//
24 template<typename T1, typename T2, bool is_const_T2, typename ExprSpec>
25 class SmoothMaxOp
26 {
27 };
28 
29 // Specialization for two Sacado expressions
30 template<typename T1, typename T2>
31 class SmoothMaxOp<T1, T2, false, Sacado::Fad::Exp::ExprSpecDefault>
32  : public Expr<SmoothMaxOp<T1, T2, false, Sacado::Fad::Exp::ExprSpecDefault>>
33 {
34  public:
35  using ExprT1 = typename std::remove_cv<T1>::type;
36  using ExprT2 = typename std::remove_cv<T2>::type;
37  using value_type1 = typename ExprT1::value_type;
38  using value_type2 = typename ExprT2::value_type;
39  using value_type = typename Sacado::Promote<value_type1, value_type2>::type;
40  using scalar_type1 = typename ExprT1::scalar_type;
41  using scalar_type2 = typename ExprT2::scalar_type;
42  using scalar_type =
43  typename Sacado::Promote<scalar_type1, scalar_type2>::type;
44  using expr_spec_type = Sacado::Fad::Exp::ExprSpecDefault;
45 
46  SACADO_INLINE_FUNCTION explicit SmoothMaxOp(const T1& x,
47  const T2& y,
48  double tol)
49  : x_(x)
50  , y_(y)
51  , tol_(tol)
52  {
53  }
54 
55  SACADO_INLINE_FUNCTION int size() const
56  {
57  const int sz1 = x_.size();
58  const int sz2 = y_.size();
59  return sz1 > sz2 ? sz1 : sz2;
60  }
61 
62  SACADO_INLINE_FUNCTION bool hasFastAccess() const
63  {
64  return x_.hasFastAccess() && y_.hasFastAccess();
65  }
66 
67  SACADO_INLINE_FUNCTION value_type val() const
68  {
69  if (x_ <= y_ - tol_)
70  {
71  return y_.val();
72  }
73  else if (x_ >= y_ + tol_)
74  {
75  return x_.val();
76  }
77  else
78  {
79  return 0.5
80  * (x_.val() + y_.val()
81  + 0.5
82  * ((x_.val() - y_.val()) * (x_.val() - y_.val())
83  / tol_
84  + tol_));
85  }
86  }
87 
88  SACADO_INLINE_FUNCTION value_type dx(int i) const
89  {
90  if (x_ <= y_ - tol_)
91  {
92  return y_.dx(i);
93  }
94  else if (x_ >= y_ + tol_)
95  {
96  return x_.dx(i);
97  }
98  else
99  {
100  return 0.5
101  * (x_.dx(i) + y_.dx(i)
102  + (x_.val() - y_.val()) * (x_.dx(i) - y_.dx(i)) / tol_);
103  }
104  }
105 
106  SACADO_INLINE_FUNCTION value_type fastAccessDx(int i) const
107  {
108  if (x_ <= y_ - tol_)
109  {
110  return y_.fastAccessDx(i);
111  }
112  else if (x_ >= y_ + tol_)
113  {
114  return x_.fastAccessDx(i);
115  }
116  else
117  {
118  return 0.5
119  * (x_.fastAccessDx(i) + y_.fastAccessDx(i)
120  + (x_.val() - y_.val())
121  * (x_.fastAccessDx(i) - y_.fastAccessDx(i)) / tol_);
122  }
123  }
124 
125  private:
126  const T1& x_;
127  const T2& y_;
128  double tol_;
129 };
130 
131 // Specialization for Sacado expression, constant
132 template<typename T1, typename T2>
133 class SmoothMaxOp<T1, T2, true, Sacado::Fad::Exp::ExprSpecDefault>
134  : public Expr<SmoothMaxOp<T1, T2, true, Sacado::Fad::Exp::ExprSpecDefault>>
135 {
136  public:
137  using ExprT1 = typename std::remove_cv<T1>::type;
138  using ConstT = T2;
139  using value_type = typename ExprT1::value_type;
140  using scalar_type = typename ExprT1::scalar_type;
141  using expr_spec_type = Sacado::Fad::Exp::ExprSpecDefault;
142 
143  SACADO_INLINE_FUNCTION explicit SmoothMaxOp(const T1& x,
144  const ConstT& c,
145  double tol)
146  : x_(x)
147  , c_(c)
148  , tol_(tol)
149  {
150  }
151 
152  SACADO_INLINE_FUNCTION int size() const { return x_.size(); }
153 
154  SACADO_INLINE_FUNCTION bool hasFastAccess() const
155  {
156  return x_.hasFastAccess();
157  }
158 
159  SACADO_INLINE_FUNCTION value_type val() const
160  {
161  if (x_ <= c_ - tol_)
162  {
163  return c_;
164  }
165  else if (x_ >= c_ + tol_)
166  {
167  return x_.val();
168  }
169  else
170  {
171  return 0.5
172  * (x_.val()
173  + 0.5 * ((x_.val() - c_) * (x_.val() - c_) / tol_ + tol_));
174  }
175  }
176 
177  SACADO_INLINE_FUNCTION value_type dx(int i) const
178  {
179  if (x_ <= c_ - tol_)
180  {
181  return 0.0;
182  }
183  else if (x_ >= c_ + tol_)
184  {
185  return x_.dx(i);
186  }
187  else
188  {
189  return 0.5 * (x_.dx(i) + (x_.val() - c_) * x_.dx(i) / tol_);
190  }
191  }
192 
193  SACADO_INLINE_FUNCTION value_type fastAccessDx(int i) const
194  {
195  if (x_ <= c_ - tol_)
196  {
197  return 0.0;
198  }
199  else if (x_ >= c_ + tol_)
200  {
201  return x_.fastAccessDx(i);
202  }
203  else
204  {
205  return 0.5
206  * (x_.fastAccessDx(i)
207  + (x_.val() - c_) * x_.fastAccessDx(i) / tol_);
208  }
209  }
210 
211  private:
212  const T1& x_;
213  const ConstT& c_;
214  double tol_;
215 };
216 
217 // Generic implementation for POD types
218 template<typename T>
219 SACADO_INLINE_FUNCTION
220  typename std::enable_if<std::is_trivial<T>::value, T>::type
221  max(const T x, const T y, const double tol)
222 {
223  if (x <= y - tol)
224  {
225  return y;
226  }
227  else if (x >= y + tol)
228  {
229  return x;
230  }
231  else
232  {
233  return 0.5 * (x + y + 0.5 * ((x - y) * (x - y) / tol + tol));
234  }
235 }
236 
237 using Sacado::Fad::Exp::Expr;
238 using Sacado::Fad::Exp::ExprLevel;
239 using Sacado::Fad::Exp::IsFadExpr;
240 
241 // Implementation when both arguments are FAD expressions
242 // See Sacado_SFINAE_Macros.hpp: SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR
243 template<typename T1, typename T2>
244 SACADO_INLINE_FUNCTION typename Sacado::mpl::enable_if_c<
245  IsFadExpr<T1>::value && IsFadExpr<T2>::value
246  && ExprLevel<T1>::value == ExprLevel<T2>::value,
247  SmoothMaxOp<typename Expr<T1>::derived_type,
248  typename Expr<T2>::derived_type,
249  false,
250  typename T1::expr_spec_type>>::type
251 max(const T1& x, const T2& y, const double tol)
252 {
253  using expr_t = SmoothMaxOp<typename Expr<T1>::derived_type,
254  typename Expr<T2>::derived_type,
255  false,
256  typename T1::expr_spec_type>;
257  return expr_t(x.derived(), y.derived(), tol);
258 }
259 
260 // Specialization when first argument is constant, second is FAD,
261 // and the FAD value_type and scalar_type match
262 template<typename T>
263 SACADO_INLINE_FUNCTION SmoothMaxOp<typename Expr<T>::derived_type,
264  typename T::value_type,
265  true,
266  typename T::expr_spec_type>
267 max(const typename T::value_type& c, const Expr<T>& y, const double tol)
268 {
269  using ConstT = typename T::value_type;
270  using expr_t = SmoothMaxOp<typename Expr<T>::derived_type,
271  ConstT,
272  true,
273  typename T::expr_spec_type>;
274 
275  // Reverse the order of the arguments
276  return expr_t(y.derived(), c, tol);
277 }
278 
279 // Specialization when first argument is FAD, second is constant,
280 // and the FAD value_type and scalar_type match
281 template<typename T>
282 SACADO_INLINE_FUNCTION SmoothMaxOp<typename Expr<T>::derived_type,
283  typename T::value_type,
284  true,
285  typename T::expr_spec_type>
286 max(const Expr<T>& x, const typename T::value_type& c, const double tol)
287 {
288  using ConstT = typename T::value_type;
289  using expr_t = SmoothMaxOp<typename Expr<T>::derived_type,
290  ConstT,
291  true,
292  typename T::expr_spec_type>;
293  return expr_t(x.derived(), c, tol);
294 }
295 
296 // Specialization when first argument is constant, second is FAD,
297 // and the FAD value_type and scalar_type do NOT match
298 // See Sacado_SFINAE_Macros.hpp: SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR
299 template<typename T>
300 SACADO_INLINE_FUNCTION typename Sacado::mpl::disable_if<
301  std::is_same<typename T::value_type, typename T::scalar_type>,
302  SmoothMaxOp<typename Expr<T>::derived_type,
303  typename T::scalar_type,
304  true,
305  typename T::expr_spec_type>>::type
306 max(const typename T::scalar_type& c, const Expr<T>& y, const double tol)
307 {
308  using ConstT = typename T::scalar_type;
309  using expr_t = SmoothMaxOp<typename Expr<T>::derived_type,
310  ConstT,
311  true,
312  typename T::expr_spec_type>;
313 
314  // Reverse the order of the arguments
315  return expr_t(y.derived(), c, tol);
316 }
317 
318 // Specialization when first argument is FAD, second is constant,
319 // and the FAD value_type and scalar_type do NOT match
320 // See Sacado_SFINAE_Macros.hpp: SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR
321 template<typename T>
322 SACADO_INLINE_FUNCTION typename Sacado::mpl::disable_if<
323  std::is_same<typename T::value_type, typename T::scalar_type>,
324  SmoothMaxOp<typename Expr<T>::derived_type,
325  typename T::scalar_type,
326  true,
327  typename T::expr_spec_type>>::type
328 max(const Expr<T>& x, const typename T::scalar_type& c, const double tol)
329 {
330  using ConstT = typename T::scalar_type;
331  using expr_t = SmoothMaxOp<typename Expr<T>::derived_type,
332  ConstT,
333  true,
334  typename T::expr_spec_type>;
335  return expr_t(x.derived(), c, tol);
336 }
337 
338 } // namespace SmoothMath
339 } // namespace VertexCFD
340 
341 namespace Sacado
342 {
343 
344 using VertexCFD::SmoothMath::SmoothMaxOp;
345 
346 namespace Fad
347 {
348 namespace Exp
349 {
350 
351 template<typename T1, typename T2, bool c2, typename E>
352 struct ExprLevel<SmoothMaxOp<T1, T2, c2, E>>
353 {
354  static constexpr unsigned value_1 = ExprLevel<T1>::value;
355  static constexpr unsigned value_2 = ExprLevel<T2>::value;
356  static constexpr unsigned value = value_1 >= value_2 ? value_1 : value_2;
357 };
358 
359 template<typename T1, typename T2, bool c2, typename E>
360 struct IsFadExpr<SmoothMaxOp<T1, T2, c2, E>>
361 {
362  static const unsigned value = true;
363 };
364 
365 } // namespace Exp
366 } // namespace Fad
367 
368 template<typename T1, typename T2, bool c2, typename E>
369 struct IsExpr<VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>>
370 {
371  static const unsigned value = true;
372 };
373 
374 template<typename T1, typename T2, bool c2, typename E>
375 struct BaseExprType<VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>>
376 {
377  using base_expr_1 = typename BaseExprType<T1>::type;
378  using base_expr_2 = typename BaseExprType<T2>::type;
379  using type = typename Promote<base_expr_1, base_expr_2>::type;
380 };
381 
382 template<typename T1, typename T2, bool c2, typename E>
383 struct IsSimdType<VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>>
384 {
385  static const bool value
386  = IsSimdType<typename VertexCFD::SmoothMath::
387  SmoothMaxOp<T1, T2, c2, E>::scalar_type>::value;
388 };
389 
390 template<typename T1, typename T2, bool c2, typename E>
391 struct ValueType<VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>>
392 {
393  using type =
394  typename VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>::value_type;
395 };
396 
397 template<typename T1, typename T2, bool c2, typename E>
398 struct ScalarType<VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>>
399 {
400  using type =
401  typename VertexCFD::SmoothMath::SmoothMaxOp<T1, T2, c2, E>::scalar_type;
402 };
403 
404 } // namespace Sacado
405 
406 #endif // end VERTEXCFD_UTILS_SMOOTHMAX_HPP
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
VertexCFD::SmoothMath::SmoothMaxOp
Definition: VertexCFD_Utils_SmoothMax.hpp:26
Generated by doxygen 1.8.20