OpenLB 1.7
Loading...
Searching...
No Matches
aDiff.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2012-2021 Mathias J. Krause, Jan E. Marquardt, Luca Heim
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
29#ifndef A_DIFF_H
30#define A_DIFF_H
31
32#include <cstdlib>
33#include <type_traits>
34#include <array>
35#include <iostream>
36#include <cassert>
37#include <type_traits>
38
39#include "core/baseType.h"
40#include "core/vector.h"
41#include "core/util.h"
42
43struct AD { };
44
45// All OpenLB code is contained in this namespace.
46namespace olb {
47
48namespace util {
49
52
58// Used to enable ADf instances in scalar-vector operations
63template <class T, unsigned DIM>
64class ADf final : public AD {
65private:
66 inline constexpr void checkDataType();
67public:
68 using base_t = T;
69 static constexpr unsigned dim = DIM;
70
72 T _v = T();
74 Vector<T,DIM> _d = ( T( 0 ) );
75
76 // class functions
77 inline constexpr ADf();
78
79 inline constexpr ADf(const ADf& a);
80 inline constexpr ADf(const T& v);
81 inline constexpr ADf(const T& v, const Vector<T,DIM>& d);
82
83 constexpr ADf(ADf&& a);
84
85 inline constexpr T& v();
86 inline constexpr T& d(unsigned i);
87 inline constexpr const T& d(unsigned i) const;
88 inline constexpr Vector<T,DIM>& d();
89
90 inline constexpr ADf& operator = (const ADf& a);
91 inline constexpr ADf& operator = (ADf&& a);
92 inline constexpr ADf& operator = (const T& v);
93
94 inline constexpr ADf& operator += (const ADf& a);
95 inline constexpr ADf& operator += (const T& v);
96
97 inline constexpr ADf& operator -= (const ADf& a);
98 inline constexpr ADf& operator -= (const T& v);
99
100 inline constexpr ADf& operator *= (const ADf& a);
101 inline constexpr ADf& operator *= (const T& v);
102
103 inline constexpr ADf& operator /= (const ADf& a);
104 inline constexpr ADf& operator /= (const T& v);
105
106 inline constexpr void setDiffVariable(unsigned iD);
107
108 inline constexpr bool hasZeroDerivative() const;
109 inline constexpr operator base_t() const;
110};
111
112
113template< class T >
114struct is_adf {
115 static constexpr bool value = false;
116};
117
118template< class S, unsigned DIM >
119struct is_adf<ADf<S,DIM>> {
120 static constexpr bool value = true;
121};
122
123template< class T >
124inline constexpr bool is_adf_v = is_adf<T>::value;
125
126// class functions
127template <class T, unsigned DIM>
128inline constexpr ADf<T,DIM>::ADf():ADf<T,DIM>( T() )
129{ }
130
131template <class T, unsigned DIM>
132inline constexpr ADf<T,DIM>::ADf(const ADf& a):ADf<T,DIM>(a._v, a._d)
133{ }
134
135template <class T, unsigned DIM>
136inline constexpr ADf<T,DIM>::ADf(const T& v):ADf<T,DIM>(v,Vector<T,DIM> {})
137{ }
138
139template <class T, unsigned DIM>
140inline constexpr ADf<T,DIM>::ADf(const T& v, const Vector<T,DIM>& d): _v(v), _d(d)
141{
142 checkDataType();
143}
144
145template <class T, unsigned DIM>
146inline constexpr ADf<T,DIM>::ADf(ADf&& a):ADf<T,DIM>(a._v, a._d)
147{ }
148
150template <class T, unsigned DIM>
151inline constexpr void ADf<T,DIM>::checkDataType()
152{
153 static_assert(std::is_floating_point<BaseType<T>>::value,
154 "Use floating point types (float, double or long double). Do not use integer types.");
155}
156
157template <class T, unsigned DIM>
158inline constexpr ADf<T,DIM>& ADf<T,DIM>::operator= (const ADf& a)
159{
160 _v=a._v;
161 _d=a._d;
162 return *this;
163}
164
165template <class T, unsigned DIM>
167{
168 _v=a._v;
169 _d=a._d;
170 return *this;
171}
172
173template <class T, unsigned DIM>
174inline constexpr ADf<T,DIM>& ADf<T,DIM>::operator= (const T& v)
175{
176 _v=v;
177 _d = ( T() );
178 return *this;
179}
180
181template <class T, unsigned DIM>
183{
184 _v+=a._v;
185 _d+=a._d;
186 return *this;
187}
188
189template <class T, unsigned DIM>
190inline constexpr ADf<T,DIM>& ADf<T,DIM>::operator += (const T& v)
191{
192 _v+=v;
193 return *this;
194}
195
196template <class T, unsigned DIM>
198{
199 _v-=a._v;
200 _d-=a._d;
201 return *this;
202}
203
204template <class T, unsigned DIM>
205inline constexpr ADf<T,DIM>& ADf<T,DIM>::operator -= (const T& v)
206{
207 _v-=v;
208 return *this;
209}
210
211template <class T, unsigned DIM>
213{
214 _d=_d*a._v+_v*a._d;
215 _v*=a._v;
216 return *this;
217}
218
219template <class T, unsigned DIM>
220inline constexpr ADf<T,DIM>& ADf<T,DIM>::operator *= (const T& v)
221{
222 _v*=v;
223 _d*=v;
224 return *this;
225}
226
227template <class T, unsigned DIM>
229{
230 T tmp(T(1)/a._v);
231 _v*=tmp;
232 _d=tmp*(_d-_v*a._d);
233 return *this;
234}
235
236template <class T, unsigned DIM>
237inline constexpr ADf<T,DIM>& ADf<T,DIM>::operator /= (const T& v)
238{
239 T tmp(T(1)/v);
240 _v*=tmp;
241 _d*=tmp;
242 return *this;
243}
244
245template <class T, unsigned DIM>
246inline constexpr T& ADf<T,DIM>::v()
247{
248 return _v;
249}
250
251template <class T, unsigned DIM>
252inline constexpr T& ADf<T,DIM>::d(unsigned i)
253{
254 return _d[i];
255}
256
257template <class T, unsigned DIM>
258inline constexpr const T& ADf<T,DIM>::d(unsigned i) const
259{
260 return _d[i];
261}
262
263template <class T, unsigned DIM>
264inline constexpr Vector<T,DIM>& ADf<T,DIM>::d()
265{
266 return _d;
267}
268
269template <class T, unsigned DIM>
270inline constexpr void ADf<T,DIM>::setDiffVariable(unsigned iD)
271{
272 _d =( T(0) );
273 _d[iD] = T(1);
274}
275
276template <class T, unsigned DIM>
277inline std::ostream& operator << (std::ostream& os, const ADf<T,DIM>& o)
278{
279 os << o._v;
280 if constexpr (DIM>0) {
281 os << o._d;
282 }
283 return os;
284}
285
286// WARNING: Reads only the value, not the derivatives
287template <class T, unsigned DIM>
288inline std::istream& operator >> (std::istream& is, ADf<T,DIM>& in)
289{
290 is >> in._v;
291 in._d = T{0};
292 return is;
293}
294
295
296// indirect helper functions, providing general arithmetics
297// overloading of operators for interaction with BaseType and itself
298// addition, subtraction, multiplication and division
299
300template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
301inline constexpr ADf<T,DIM> operator+ (const U& a, const ADf<T,DIM>& b)
302{
303 return ADf<T,DIM>(b) += olb::BaseType<ADf<T,DIM>>(a);
304}
305
306template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
307inline constexpr ADf<T,DIM> operator+ (const ADf<T,DIM>& a,const U& b)
308{
309 return ADf<T,DIM>(a) += olb::BaseType<ADf<T,DIM>>(b);
310}
311
312template <class T, unsigned DIM>
313inline constexpr ADf<T,DIM> operator+ (const ADf<T,DIM>& a, const ADf<T,DIM>& b)
314{
315 return ADf<T,DIM>(a) += b;
316}
317
318template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
319inline constexpr ADf<T,DIM> operator- (const U& a, const ADf<T,DIM>& b)
320{
322}
323
324template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
325inline constexpr ADf<T,DIM> operator- (const ADf<T,DIM>& a,const U& b)
326{
327 return ADf<T,DIM>(a) -= olb::BaseType<ADf<T,DIM>>(b);
328}
329
330template <class T, unsigned DIM>
331inline constexpr ADf<T,DIM> operator- (const ADf<T,DIM>& a, const ADf<T,DIM>& b)
332{
333 return ADf<T,DIM>(a)-=b;
334}
335
336template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
337inline constexpr ADf<T,DIM> operator* (const U& a, const ADf<T,DIM>& b)
338{
339 return ADf<T,DIM>(b) *= olb::BaseType<ADf<T,DIM>>(a);
340}
341
342template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
343inline constexpr ADf<T,DIM> operator* (const ADf<T,DIM>& a,const U& b)
344{
345 return ADf<T,DIM>(a) *= olb::BaseType<ADf<T,DIM>>(b);
346}
347
348template <class T, unsigned DIM>
349inline constexpr ADf<T,DIM> operator* (const ADf<T,DIM>& a, const ADf<T,DIM>& b)
350{
351 return ADf<T,DIM>(a) *= b;
352}
353
354template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
355inline constexpr ADf<T,DIM> operator/ (const U& a, const ADf<T,DIM>& b)
356{
357 return ADf<T,DIM>(a) /= b;
358}
359
360template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
361inline constexpr ADf<T,DIM> operator/ (const ADf<T,DIM>& a,const U& b)
362{
363 return ADf<T,DIM>(a) /= olb::BaseType<ADf<T,DIM>>(b);
364}
365
366
367template <class T, unsigned DIM>
368inline constexpr ADf<T,DIM> operator/ (const ADf<T,DIM>& a, const ADf<T,DIM>& b)
369{
370 return ADf<T,DIM>(a) /= b;
371}
372
373/* Unary operators */
374template <class T, unsigned DIM>
375inline constexpr ADf<T,DIM> operator+ (const ADf<T,DIM>& a)
376{
377 return ADf<T,DIM>(a);
378}
379
380template <class T, unsigned DIM>
381inline constexpr ADf<T,DIM> operator- (const ADf<T,DIM>& a)
382{
383 return ADf<T,DIM>(a) *= -1;
384}
385
386
387
388
389
390
391
392
393// extended math functions:
394// pow, sqr, exp, log, log10, sqrt, sin, cos, tan, asin, acos, atan
395template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
396inline ADf<T,DIM> pow (const U& a, const ADf<T,DIM>& b)
397{
398 ADf<T,DIM> c(std::pow(T(a),b._v), b._d);
399 T tmp(c._v*std::log(T(a)));
400 c._d*=tmp;
401 return c;
402}
403
404template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
405inline ADf<T,DIM> pow (const ADf<T,DIM>& a,const U& b)
406{
407 ADf<T,DIM> c(std::pow(a._v,T(b)), a._d);
408 T tmp(T(b)*std::pow(a._v,b-T(1)));
409 c._d*=tmp;
410 return c;
411}
412
413
414template <class T, unsigned DIM>
415inline ADf<T,DIM> pow (const ADf<T,DIM>& a, const ADf<T,DIM>& b)
416{
417 if (!util::nearZero(a._v)) {
418 return exp(b*log(a));
419 }
420 else {
421 ADf<T,DIM> c(std::pow(a._v, b._v), Vector<T,DIM> {std::numeric_limits<T>::quiet_NaN()});
422#ifdef AdWarnings
423 std::cout << "ADf WARNING: pow(Adf) - pow evaluated at non-differentiable point" << std::endl;
424#endif
425 if (std::isfinite(c._v)) {
426 for (unsigned i = 0; i < DIM; ++i) {
427 if (util::nearZero(a.d(i))) {
428 c._d[i] = T(0);
429 }
430 }
431 }
432 return c;
433 }
434}
435
436template <class T, unsigned DIM>
437inline ADf<T,DIM> pow (const ADf<T,DIM>& a, int b)
438{
439 ADf<T,DIM> c(std::pow(a._v,b), a._d);
440 T tmp(b*std::pow(a._v, b-1));
441 c._d*=tmp;
442 return c;
443}
444
445template <class T, unsigned DIM>
446inline ADf<T,DIM> sqr (const ADf<T,DIM>& a)
447{
448 ADf<T,DIM> c((a._v)*(a._v), a._d);
449 T tmp(T(2)*a._v);
450 c._d*=tmp;
451 return c;
452}
453
454template <class T, unsigned DIM>
455inline ADf<T,DIM> exp (const ADf<T,DIM>& a)
456{
457 ADf<T,DIM> c(std::exp(a._v), a._d);
458 c._d*=c._v;
459 return c;
460}
461
462template <unsigned DIM>
464{
465 return exp(a);
466}
467
468template <unsigned DIM>
470{
471 return exp(a);
472}
473
474template <class T, unsigned DIM>
475inline ADf<T,DIM> log (const ADf<T,DIM>& a)
476{
477 ADf<T,DIM> c(std::log(a._v), a._d);
478 T tmp (T(1)/a._v);
479 c._d*=tmp;
480 return c;
481}
482
483template <unsigned DIM>
485{
486 return log(a);
487}
488
489template <unsigned DIM>
491{
492 return log(a);
493}
494
495template <class T, unsigned DIM>
496inline ADf<T,DIM> log10 (const ADf<T,DIM>& a)
497{
498 ADf<T,DIM> c(std::log10(a._v), a._d);
499 T tmp (T(1)/(a._v*std::log(T(10))));
500 c._d*=tmp;
501 return c;
502}
503
504template <unsigned DIM>
506{
507 return log10(a);
508}
509
510template <unsigned DIM>
512{
513 return log10(a);
514}
515
516template <class T, unsigned DIM>
517inline ADf<T,DIM> log2 (const ADf<T,DIM>& a)
518{
519 ADf<T,DIM> c(std::log2(a._v), a._d);
520 T tmp (T(1)/(a._v*std::log(T(2))));
521 c._d*=tmp;
522 return c;
523}
524
525template <unsigned DIM>
527{
528 return log2(a);
529}
530
531template <unsigned DIM>
533{
534 return log2(a);
535}
536
537template <class T, unsigned DIM>
538inline ADf<T,DIM> log1p (const ADf<T,DIM>& a)
539{
540 return log(a+T(1));
541}
542
543template <unsigned DIM>
545{
546 return log1p(a);
547}
548
549template <unsigned DIM>
551{
552 return log1p(a);
553}
554
555template <class T, unsigned DIM>
556inline ADf<T,DIM> sqrt (const ADf<T,DIM>& a)
557{
558 ADf<T,DIM> c(std::sqrt(a._v), a._d);
559 T tmp(T(1.)/(c._v*T(2)));
560 for (unsigned i = 0; i < DIM; ++i) {
561 if (!util::nearZero(a.d(i))) {
562 c._d[i] *= tmp;
563 }
564 }
565 return c;
566}
567
568template <class T, unsigned DIM>
569inline ADf<T,DIM> sin (const ADf<T,DIM>& a)
570{
571 ADf<T,DIM> c(std::sin(a._v), a._d);
572 T tmp(std::cos(a._v));
573 c._d*=tmp;
574 return c;
575}
576
577template <class T, unsigned DIM>
578inline ADf<T,DIM> cos (const ADf<T,DIM>& a)
579{
580 ADf<T,DIM> c(std::cos(a._v), a._d);
581 T tmp(-std::sin(a._v));
582 c._d*=tmp;
583 return c;
584}
585
586template <class T, unsigned DIM>
587inline ADf<T,DIM> tan (const ADf<T,DIM>& a)
588{
589 ADf<T,DIM> c(std::tan(a._v), a._d);
590 T tmp(T(1)+c._v*c._v);
591 c._d*=tmp;
592 return c;
593}
594
595template <class T, unsigned DIM>
596inline ADf<T,DIM> asin (const ADf<T,DIM>& a)
597{
598 ADf<T,DIM> c(std::asin(a._v), a._d);
599 T tmp(T(1)/std::sqrt(T(1)-(a._v)*(a._v)));
600 c._d*=tmp;
601 return c;
602}
603
604template <class T, unsigned DIM>
605inline ADf<T,DIM> acos (const ADf<T,DIM>& a)
606{
607 ADf<T,DIM> c(std::acos(a._v), a._d);
608 T tmp(-T(1)/std::sqrt(T(1)-(a._v)*(a._v)));
609 c._d*=tmp;
610 return c;
611}
612
613template <class T, unsigned DIM>
614inline ADf<T,DIM> atan (const ADf<T,DIM>& a)
615{
616 ADf<T,DIM> c(std::atan(a._v), a._d);
617 T tmp(T(1)/(T(1)+(a._v)*(a._v)));
618 c._d*=tmp;
619 return c;
620}
621
622template <class T, unsigned DIM>
623inline ADf<T,DIM> atan2 (const T& y, const ADf<T,DIM>& x)
624{
625 ADf<T,DIM> c(std::atan2(y, x._v));
626 T tmpB(-y / (x._v*x._v + y*y));
627 c._d = tmpB*x._d;
628 return c;
629}
630
631template <class T, unsigned DIM>
632inline ADf<T,DIM> atan2 (const ADf<T,DIM>& y, const T& x)
633{
634 ADf<T,DIM> c(std::atan2(y._v, x));
635 T tmpA(x / (x*x + y._v*y._v));
636 c._d = tmpA*y._d;
637 return c;
638}
639
640template <class T, unsigned DIM>
641inline ADf<T,DIM> atan2 (const ADf<T,DIM>& y, const ADf<T,DIM>& x)
642{
643 ADf<T,DIM> c(std::atan2(y._v, x._v));
644 T tmpA(x._v / (x._v*x._v + y._v*y._v));
645 T tmpB(-y._v / (x._v*x._v + y._v*y._v));
646 c._d = tmpA*y._d + tmpB*x._d;
647 return c;
648}
649
650template <class T, unsigned DIM>
651inline ADf<T,DIM> sinh (const ADf<T,DIM>& a)
652{
653 return 0.5*(exp(a)-exp(-a));
654}
655
656template <class T, unsigned DIM>
657inline ADf<T,DIM> cosh (const ADf<T,DIM>& a)
658{
659 return 0.5*(exp(a)+exp(-a));
660}
661
662template <class T, unsigned DIM>
663inline ADf<T,DIM> tanh (const ADf<T,DIM>& a)
664{
665 return 1 - 2.0 / (exp(2*a) + 1);
666}
667
668template <class T, unsigned DIM>
669inline ADf<T,DIM> asinh (const ADf<T,DIM>& a)
670{
671 return log( a + sqrt(a*a+1) );
672}
673
674template <class T, unsigned DIM>
675inline ADf<T,DIM> acosh (const ADf<T,DIM>& a)
676{
677 if (a._v >= 1) {
678 return log( a + sqrt(a*a-1) );
679 }
680 else {
681 ADf<T,DIM> c;
682 c._v = std::acosh(a._v);
683 c._d = ( std::numeric_limits<T>::quiet_NaN() );
684 return c;
685 }
686}
687
688template <class T, unsigned DIM>
689inline ADf<T,DIM> atanh (const ADf<T,DIM>& a)
690{
691 if (std::abs(a._v) < 1) {
692 return 0.5 * log( (1+a) / (1-a) );
693 }
694 else {
695 ADf<T,DIM> c;
696 c._v = std::atanh(a._v);
697 c._d = ( std::numeric_limits<T>::quiet_NaN() );
698 return c;
699 }
700}
701
702template <class T, unsigned DIM>
703inline ADf<T,DIM> fmod (const ADf<T,DIM>& a, const ADf<T,DIM>& b)
704{
705 const T k = std::floor(T(a) / T(b));
706 return ADf<T,DIM>(a - k * b);
707}
708
709template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
710inline ADf<T,DIM> fmod (const ADf<T,DIM>& a, const U& b)
711{
712 return fmod(a, ADf<T,DIM>(b));
713}
714
715template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
716inline ADf<T,DIM> fmod (const U& a, const ADf<T,DIM>& b)
717{
718 return fmod( ADf<T,DIM>(a),b);
719}
720
721
722
723
724
725
726
727
728
730template <class T, unsigned DIM>
731constexpr inline bool ADf<T,DIM>::hasZeroDerivative() const
732{
733 T sum = 0;
734 for (unsigned i=0; i<DIM; i++) {
735 sum += std::fabs(_d[i]);
736 }
737 return sum==0;
738}
739
740
741
742// boolean comparison operators: ==, !=, >, >=, <, <=
743template <class T, unsigned DIM>
744inline constexpr bool operator == (const ADf<T,DIM> &a, const ADf<T,DIM> &b)
745{
746 return a._v==b._v;
747}
748
749template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
750inline constexpr bool operator == (const ADf<T,DIM> &a, const U &b)
751{
752 return a._v==olb::BaseType<ADf<T,DIM>>(b);
753}
754
755template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
756inline constexpr bool operator == (const U &a, const ADf<T,DIM> &b)
757{
758 return olb::BaseType<ADf<T,DIM>>(a)==b._v;
759}
760
761template <class T, unsigned DIM>
762inline constexpr bool operator != (const ADf<T,DIM> &a, const ADf<T,DIM> &b)
763{
764 return a._v!=b._v;
765}
766
767template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
768inline constexpr bool operator != (const ADf<T,DIM> &a, const U &b)
769{
770 return a._v!=olb::BaseType<ADf<T,DIM>>(b);
771}
772
773template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
774inline constexpr bool operator != (const U &a, const ADf<T,DIM> &b)
775{
776 return olb::BaseType<ADf<T,DIM>>(a)!=b._v;
777}
778
779
780template <class T, unsigned DIM>
781inline constexpr bool operator > (const ADf<T,DIM> &a, const ADf<T,DIM> &b)
782{
783 return a._v>b._v;
784}
785
786template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
787inline constexpr bool operator > (const ADf<T,DIM> &a, const U &b)
788{
789 return a._v>olb::BaseType<ADf<T,DIM>>(b);
790}
791
792template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
793inline constexpr bool operator > (const U &a, const ADf<T,DIM> &b)
794{
795 return olb::BaseType<ADf<T,DIM>>(a)>b._v;
796}
797
798template <class T, unsigned DIM>
799inline constexpr bool operator >= (const ADf<T,DIM> &a, const ADf<T,DIM> &b)
800{
801 return a._v>=b._v;
802}
803
804template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
805inline constexpr bool operator >= (const ADf<T,DIM> &a, const U &b)
806{
807 return a._v>=olb::BaseType<ADf<T,DIM>>(b);
808}
809
810template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
811inline constexpr bool operator >= (const U &a, const ADf<T,DIM> &b)
812{
813 return olb::BaseType<ADf<T,DIM>>(a)>=b._v;
814}
815
816
817template <class T, unsigned DIM>
818inline constexpr bool operator < (const ADf<T,DIM> &a, const ADf<T,DIM> &b)
819{
820 return a._v<b._v;
821}
822
823template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
824inline constexpr bool operator < (const ADf<T,DIM> &a, const U &b)
825{
826 return a._v<olb::BaseType<ADf<T,DIM>>(b);
827}
828
829template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
830inline constexpr bool operator < (const U &a, const ADf<T,DIM> &b)
831{
832 return olb::BaseType<ADf<T,DIM>>(a)<b._v;
833}
834
835template <class T, unsigned DIM>
836inline constexpr bool operator <= (const ADf<T,DIM> &a, const ADf<T,DIM> &b)
837{
838 return a._v<=b._v;
839}
840
841template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
842inline constexpr bool operator <= (const ADf<T,DIM> &a, const U &b)
843{
844 return a._v<=olb::BaseType<ADf<T,DIM>>(b);
845}
846
847template <class T, unsigned DIM, typename U, typename std::enable_if<std::is_integral<U>::value | std::is_floating_point<U>::value,int>::type = 0>
848inline constexpr bool operator <= (const U &a, const ADf<T,DIM> &b)
849{
850 return olb::BaseType<ADf<T,DIM>>(a)<=b._v;
851}
852
853// typecast functions
855template <class T, unsigned DIM>
856inline constexpr ADf<T,DIM>::operator ADf<T,DIM>::base_t() const
857{
858#ifdef AdWarnings
859 if (this->hasZeroDerivative() != true) {
860 std::cout << "ADf WARNING: non-zero derivative in typecast ADf -> Basetype" << std::endl;
861 std::cout << *this << std::endl;
862 //exit(1);
863 }
864#endif
865 return base_t(_v);
866}
867
868template <class T, unsigned DIM>
869inline ADf<T,DIM> floor (const ADf<T,DIM>& a)
870{
871 ADf<T,DIM> c;
872 if (a._v == std::floor(a._v)) {
873 c._v = a._v;
874#ifdef AdWarnings
875 std::cout << "ADf WARNING: floor(Adf) - floor evaluated at non-differentiable point" << std::endl;
876#endif
877 for (unsigned i = 0; i < DIM; ++i) {
878 c._d[i] = (a.d(i) == 0) ? T(0) : std::numeric_limits<T>::quiet_NaN();
879 }
880 return c;
881 }
882 c._v = std::floor(a._v);
883 c._d = ( T(0) );
884 return c;
885}
886
887template <unsigned DIM>
889{
890 return floor(a);
891}
892
893template <unsigned DIM>
895{
896 return floor(a);
897}
898
899template <class T, unsigned DIM>
900inline ADf<T,DIM> ceil (const ADf<T,DIM>& a)
901{
902 // Similar to floor
903 ADf<T,DIM> c (std::ceil(a._v));
904 if (a._v == c._v) {
905#ifdef AdWarnings
906 std::cout << "ADf WARNING: ceil(Adf) - ceil evaluated at non-differentiable point" << std::endl;
907#endif
908 for (unsigned i = 0; i < DIM; ++i) {
909 c._d[i] = (a.d(i) == 0) ? T(0) : std::numeric_limits<T>::quiet_NaN();
910 }
911 }
912 return c;
913}
914
915template <unsigned DIM>
917{
918 return ceil(a);
919}
920
921template <unsigned DIM>
923{
924 return ceil(a);
925}
926
927template <class T, unsigned DIM>
928inline ADf<T,DIM> round (const ADf<T,DIM>& a)
929{
930 ADf<T,DIM> c (std::round(a._v));
931 constexpr T EPSILON (std::numeric_limits<T>::epsilon());
932 if (std::round(a._v - EPSILON) != std::round(a._v + EPSILON)) {
933#ifdef AdWarnings
934 std::cout << "ADf WARNING: round(Adf) - round evaluated at non-differentiable point" << std::endl;
935#endif
936 for (unsigned i = 0; i < DIM; ++i) {
937 c._d[i] = (a.d(i) == 0) ? T(0) : std::numeric_limits<T>::quiet_NaN();
938 }
939 }
940 return c;
941}
942
943template <unsigned DIM>
945{
946 return round(a);
947}
948
949template <unsigned DIM>
951{
952 return round(a);
953}
954
955// lround, llround, etc. throw error, since we do not support integers with ADf
956template <class T, unsigned DIM>
957inline ADf<T,DIM> lround (const ADf<T,DIM>& a)
958{
959#ifdef AdWarnings
960 std::cout << "ADf WARNING: ADf does not support integer types. Using round() instead." << std::endl;
961#endif
962 return round(a);
963}
964
965template <unsigned DIM>
967{
968 return lround(a);
969}
970
971template <unsigned DIM>
973{
974 return lround(a);
975}
976
977template <class T, unsigned DIM>
979{
980 return lround(a);
981}
982
983template <unsigned DIM>
985{
986 return llround(a);
987}
988
989template <unsigned DIM>
991{
992 return llround(a);
993}
994
995template <class T, unsigned DIM>
996inline ADf<T,DIM> fabs (const ADf<T,DIM>& a)
997{
998#ifdef AdWarnings
999 if (a._v == T(0)) {
1000 std::cout << "ADf WARNING: fabs(Adf) - fabs evaluated at non-differentiable point" << std::endl;
1001 }
1002#endif
1003 return ADf<T,DIM>( util::sign(a._v) * a );
1004}
1005
1006template <unsigned DIM>
1008{
1009 return fabs(a);
1010}
1011
1012template <unsigned DIM>
1014{
1015 return fabs(a);
1016}
1017
1018template <class T, unsigned DIM>
1019inline ADf<T,DIM> abs (const ADf<T,DIM>& a)
1020{
1021 return fabs(a);
1022}
1023
1024template <class T, unsigned DIM>
1025inline ADf<T,DIM> labs (const ADf<T,DIM>& a)
1026{
1027#ifdef AdWarnings
1028 std::cout << "ADf WARNING: ADf is not supporting integer types. Using abs() instead." << std::endl;
1029#endif
1030 return abs(a);
1031}
1032
1033template <class T, unsigned DIM>
1034inline ADf<T,DIM> llabs (const ADf<T,DIM>& a)
1035{
1036 return labs(a);
1037}
1038
1039template <class T, unsigned DIM>
1040inline constexpr ADf<T,DIM> max (const olb::BaseType<ADf<T,DIM>>& a, const ADf<T,DIM>& b)
1041{
1042 ADf<T,DIM> c(a);
1043 c = max(c, b);
1044 return c;
1045}
1046
1047template <class T, unsigned DIM>
1048inline constexpr ADf<T,DIM> max(const ADf<T,DIM>& a,const olb::BaseType<ADf<T,DIM>>& b)
1049{
1050 ADf<T,DIM> c(b);
1051 c = max(a, c);
1052 return c;
1053}
1054
1055template <class T, unsigned DIM>
1056inline constexpr ADf<T,DIM> max(const ADf<T,DIM>& a, const ADf<T,DIM>& b)
1057{
1058 ADf<T,DIM> c((a + b + fabs(a - b) ) * T(0.5));
1059 return c;
1060}
1061
1062template <class T, unsigned DIM>
1063inline constexpr ADf<T,DIM> min (const olb::BaseType<ADf<T,DIM>>& a, const ADf<T,DIM>& b)
1064{
1065 ADf<T,DIM> c(a);
1066 c = min(c, b);
1067 return c;
1068}
1069
1070template <class T, unsigned DIM>
1071inline constexpr ADf<T,DIM> min(const ADf<T,DIM>& a,const olb::BaseType<ADf<T,DIM>>& b)
1072{
1073 ADf<T,DIM> c(b);
1074 c = min(a, c);
1075 return c;
1076}
1077
1078template <class T, unsigned DIM>
1079inline constexpr ADf<T,DIM> min(const ADf<T,DIM>& a, const ADf<T,DIM>& b)
1080{
1081 ADf<T,DIM> c((a + b - fabs(a - b) ) * T(0.5));
1082 return c;
1083}
1084
1085// overload nearZero to prevent bad behaviour due to implicit casting
1086template <class T, unsigned DIM>
1087inline bool nearZero(const ADf<T,DIM>& a)
1088{
1089 const T EPSILON = std::numeric_limits<T>::epsilon();
1090 if (a._v > -EPSILON && a._v < EPSILON) {
1091 return true;
1092 }
1093 else {
1094 return false;
1095 }
1096}
1097
1098template <typename T>
1099inline constexpr void throwADfException(T arg)
1100{
1101 static_assert(std::is_floating_point<T>::value,
1102 "The data type ADf has slipped in here. If you really want to use it here and lose all derivation information, then perform an explicit typecast.");
1103}
1104
1105} // namespace util
1106
1107} // namespace olb
1108
1109#endif
Plain old scalar vector.
Definition vector.h:47
Definition of a description of a algoritmic differentiation data type using the forward method.
Definition aDiff.h:64
constexpr T & v()
Definition aDiff.h:246
static constexpr unsigned dim
Definition aDiff.h:69
constexpr ADf & operator/=(const ADf &a)
Definition aDiff.h:228
constexpr ADf()
Definition aDiff.h:128
constexpr void setDiffVariable(unsigned iD)
Definition aDiff.h:270
T _v
value
Definition aDiff.h:72
constexpr T & d(unsigned i)
Definition aDiff.h:252
Vector< T, DIM > _d
derivatives
Definition aDiff.h:74
constexpr bool hasZeroDerivative() const
tests if ADf has only zero derivatives
Definition aDiff.h:731
constexpr Vector< T, DIM > & d()
Definition aDiff.h:264
constexpr ADf & operator=(const ADf &a)
Definition aDiff.h:158
constexpr ADf & operator-=(const ADf &a)
Definition aDiff.h:197
constexpr ADf & operator*=(const ADf &a)
Definition aDiff.h:212
constexpr ADf & operator+=(const ADf &a)
Definition aDiff.h:182
ADf< T, DIM > labs(const ADf< T, DIM > &a)
Definition aDiff.h:1025
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
ADf< long double, DIM > fabsl(const ADf< long double, DIM > &a)
Definition aDiff.h:1013
ADf< T, DIM > abs(const ADf< T, DIM > &a)
Definition aDiff.h:1019
ADf< T, DIM > ceil(const ADf< T, DIM > &a)
Definition aDiff.h:900
ADf< T, DIM > floor(const ADf< T, DIM > &a)
Definition aDiff.h:869
ADf< T, DIM > tan(const ADf< T, DIM > &a)
Definition aDiff.h:587
ADf< T, DIM > fmod(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:703
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
ADf< T, DIM > atan2(const T &y, const ADf< T, DIM > &x)
Definition aDiff.h:623
ADf< float, DIM > llroundf(const ADf< float, DIM > &a)
Definition aDiff.h:984
ADf< T, DIM > asinh(const ADf< T, DIM > &a)
Definition aDiff.h:669
ADf< long double, DIM > log1pl(const ADf< long double, DIM > &a)
Definition aDiff.h:550
ADf< T, DIM > lround(const ADf< T, DIM > &a)
Definition aDiff.h:957
ADf< float, DIM > log10f(const ADf< float, DIM > &a)
Definition aDiff.h:505
ADf< T, DIM > log(const ADf< T, DIM > &a)
Definition aDiff.h:475
ADf< long double, DIM > expl(const ADf< long double, DIM > &a)
Definition aDiff.h:469
ADf< float, DIM > logf(const ADf< float, DIM > &a)
Definition aDiff.h:484
ADf< long double, DIM > log2l(const ADf< long double, DIM > &a)
Definition aDiff.h:532
std::ostream & operator<<(std::ostream &os, const ADf< T, DIM > &o)
Definition aDiff.h:277
constexpr ADf< T, DIM > operator*(const U &a, const ADf< T, DIM > &b)
Definition aDiff.h:337
ADf< float, DIM > roundf(const ADf< float, DIM > &a)
Definition aDiff.h:944
constexpr bool operator<=(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:836
constexpr ADf< T, DIM > operator+(const U &a, const ADf< T, DIM > &b)
Definition aDiff.h:301
ADf< T, DIM > round(const ADf< T, DIM > &a)
Definition aDiff.h:928
ADf< float, DIM > log2f(const ADf< float, DIM > &a)
Definition aDiff.h:526
ADf< long double, DIM > log10l(const ADf< long double, DIM > &a)
Definition aDiff.h:511
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
ADf< T, DIM > llround(const ADf< T, DIM > &a)
Definition aDiff.h:978
ADf< float, DIM > floorf(const ADf< float, DIM > &a)
Definition aDiff.h:888
ADf< T, DIM > sin(const ADf< T, DIM > &a)
Definition aDiff.h:569
ADf< float, DIM > expf(const ADf< float, DIM > &a)
Definition aDiff.h:463
ADf< T, DIM > acosh(const ADf< T, DIM > &a)
Definition aDiff.h:675
ADf< long double, DIM > floorl(const ADf< long double, DIM > &a)
Definition aDiff.h:894
constexpr bool operator<(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:818
ADf< float, DIM > fabsf(const ADf< float, DIM > &a)
Definition aDiff.h:1007
constexpr void throwADfException(T arg)
Definition aDiff.h:1099
ADf< T, DIM > asin(const ADf< T, DIM > &a)
Definition aDiff.h:596
ADf< T, DIM > log1p(const ADf< T, DIM > &a)
Definition aDiff.h:538
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
constexpr bool is_adf_v
Definition aDiff.h:124
ADf< float, DIM > lroundf(const ADf< float, DIM > &a)
Definition aDiff.h:966
ADf< T, DIM > cosh(const ADf< T, DIM > &a)
Definition aDiff.h:657
ADf< T, DIM > llabs(const ADf< T, DIM > &a)
Definition aDiff.h:1034
constexpr bool operator!=(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:762
ADf< T, DIM > sinh(const ADf< T, DIM > &a)
Definition aDiff.h:651
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
constexpr ADf< T, DIM > operator/(const U &a, const ADf< T, DIM > &b)
Definition aDiff.h:355
ADf< long double, DIM > lroundl(const ADf< long double, DIM > &a)
Definition aDiff.h:972
ADf< long double, DIM > roundl(const ADf< long double, DIM > &a)
Definition aDiff.h:950
ADf< T, DIM > acos(const ADf< T, DIM > &a)
Definition aDiff.h:605
ADf< long double, DIM > llroundl(const ADf< long double, DIM > &a)
Definition aDiff.h:990
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
T sqr(T arg)
Definition util.h:136
int sign(T val)
Definition util.h:54
ADf< long double, DIM > ceill(const ADf< long double, DIM > &a)
Definition aDiff.h:922
ADf< T, DIM > tanh(const ADf< T, DIM > &a)
Definition aDiff.h:663
ADf< float, DIM > log1pf(const ADf< float, DIM > &a)
Definition aDiff.h:544
ADf< T, DIM > log10(const ADf< T, DIM > &a)
Definition aDiff.h:496
constexpr bool operator>(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:781
std::istream & operator>>(std::istream &is, ADf< T, DIM > &in)
Definition aDiff.h:288
ADf< T, DIM > atan(const ADf< T, DIM > &a)
Definition aDiff.h:614
constexpr bool operator==(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:744
constexpr bool operator>=(const ADf< T, DIM > &a, const ADf< T, DIM > &b)
Definition aDiff.h:799
ADf< long double, DIM > logl(const ADf< long double, DIM > &a)
Definition aDiff.h:490
ADf< T, DIM > exp(const ADf< T, DIM > &a)
Definition aDiff.h:455
ADf< float, DIM > ceilf(const ADf< float, DIM > &a)
Definition aDiff.h:916
ADf< T, DIM > atanh(const ADf< T, DIM > &a)
Definition aDiff.h:689
constexpr ADf< T, DIM > operator-(const U &a, const ADf< T, DIM > &b)
Definition aDiff.h:319
ADf< T, DIM > log2(const ADf< T, DIM > &a)
Definition aDiff.h:517
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Definition aDiff.h:578
Top level namespace for all of OpenLB.
typename util::BaseTypeHelper< T >::type BaseType
Definition baseType.h:59
Definition aDiff.h:43
static constexpr bool value
Definition aDiff.h:115
Set of functions commonly used in LB computations – header file.
efficient implementation of a vector class