adHelpers.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2021 Julius Jeßberger
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/
 
#ifndef AD_HELPERS_H
#define AD_HELPERS_H
 
#include "aDiff.h"
 
 
namespace olb {
 
namespace util {
 
template <typename SAD>
void iniDiagonal(SAD* a, int dim) {
  for (int i = 0; i < dim; ++i) {
    a[i].setDiffVariable(i);
  }
}
 
template <typename C>
void iniDiagonal(C& a) {
  for (unsigned i = 0; i < a.size(); ++i) {
    a[i].setDiffVariable(i);
  }
}
 
 
template <typename T, typename TAD>
void copyDerivatives(T* target, const TAD* source, int length) {
  for (int j = 0; j < length; ++j) {
    for (unsigned i = 0; i < TAD::dim; ++i) {
      target[j*(TAD::dim)+i] = source[j].d(i);
    }
  }
}
 
 
template<typename T, typename S, template<typename> typename C>
C<T> copyAs(const C<S>& input) {
  C<T> result = ContainerCreator<C<T>>::create(input.size());
 
  for(std::size_t it = 0; it < input.size(); ++it) {
    result[it] = T(input[it]);
  }
  return result;
}
 
 
template<typename S>
ADf<S,1> iniAD(const S source) {
  ADf<S,1> result(source);
  result.setDiffVariable(0);
  return result;
}
 
template<unsigned n, typename S>
ADf<S,n>* iniAD(const S* source) {
  ADf<S,n>* result = new ADf<S,n>[n];
  util::copyN (result, source, n);
  util::iniDiagonal (result, n);
  return result;
}
 
template<unsigned n, typename S, template<typename> typename C>
C<ADf<S,n>> iniAD(const C<S>& source) {
  auto result = copyAs<ADf<S,n>,S,C> (source);
  util::iniDiagonal (result);
  return result;
}
 
 
 
template <typename S, typename F>
auto derivativeFAD (F f, const S input) {
  ADf<S,1> inputAD = input;
  inputAD.setDiffVariable(0);
 
  return f(inputAD).d(0);
}
 
template <unsigned sourceDIM, typename S, typename F>
auto derivativeFAD (F f, const S* input) {
  using SAD = ADf<S,sourceDIM>;
  using T = decltype(f(input));
  using TAD = ADf<T,sourceDIM>;
 
  SAD inputAD[sourceDIM];
  util::copyN(inputAD, input, sourceDIM);
  util::iniDiagonal(inputAD, sourceDIM);
 
  TAD resultAD = f(inputAD);
  T* result = new T[sourceDIM];
  copyDerivatives(result, &resultAD, 1);
  return result;
}
 
template <unsigned sourceDIM, typename T, typename S, typename F>
void derivativeFAD (F f, T* output, const S* input, unsigned targetDIM) {
  using SAD = ADf<S,sourceDIM>;
  using TAD = ADf<T,sourceDIM>;
 
  SAD inputAD[sourceDIM];
  util::copyN(inputAD, input, sourceDIM);
  util::iniDiagonal(inputAD, sourceDIM);
  TAD outputAD[targetDIM];
 
  f(outputAD, inputAD);
  copyDerivatives(output, outputAD, targetDIM);
  return;
}
 
}
 
}
 
#endif