Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/influence.hh
797f71620eb8features/netcdf-mpi

Last Change
Raw File

influence.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef INFLUENCE_HH
#define INFLUENCE_HH
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
#include <type_traits>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace influence {
namespace {
const Complex I(0, 1); /// < imaginary unit
constexpr inline Real sign(bool upper) { return (upper) ? 1. : -1.; }
}
template <bool conjugate, UInt dim_q> // dim + 1 to help dim template deduction
inline Vector<Complex, dim_q + 1> computeD(VectorProxy<const Real, dim_q>& q) {
constexpr Real s = -sign(conjugate);
const auto q_norm = q.l2norm();
return {{{Complex{0, s * q(0) / q_norm}, Complex{0, s * q(1) / q_norm},
Complex{1, 0}}}};
}
template <UInt dim_q> // dim + 1 to help dim template deduction
inline Vector<Complex, dim_q + 1> computeD2(VectorProxy<const Real, dim_q>& q) {
const auto q_norm = q.l2norm();
return {
{{Complex{0, q(0) / q_norm}, Complex{0, q(1) / q_norm}, Complex{0, 0}}}};
}
/// Class for the Kelvin tensor
template <UInt dim, UInt derivative_order>
class Kelvin;
/// See details in notebook and Marc Bonnet's notes
template <>
class Kelvin<3, 0> {
protected:
constexpr static UInt dim = 3;
public:
Kelvin(Real mu, Real nu) : mu(mu), b(4 * (1 - nu)) {}
template <bool upper, typename ST>
inline Vector<Complex, dim> applyU0(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const;
template <bool upper, typename ST>
inline Vector<Complex, dim> applyU1(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const;
protected:
const Real mu, b;
};
/* -------------------------------------------------------------------------- */
/* Base Kelvin tensor */
/* -------------------------------------------------------------------------- */
template <bool upper, typename ST> //< this is ignored
inline Vector<Complex, Kelvin<3, 0>::dim>
Kelvin<3, 0>::applyU0(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const {
auto d2 = computeD2(q);
d2 *= d2.dot(f);
d2 += b * f;
d2(2) -= f(2);
d2 *= 1. / (2. * mu * q.l2norm() * b);
return d2;
}
template <bool upper, typename ST>
inline Vector<Complex, Kelvin<3, 0>::dim>
Kelvin<3, 0>::applyU1(VectorProxy<const Real, dim - 1>& q,
StaticVector<Complex, ST, dim>& f) const {
auto d = computeD<upper>(q);
d *= d.dot(f);
d *= sign(upper) / (2. * mu * q.l2norm() * b);
return d;
}
/* -------------------------------------------------------------------------- */
/* Kelvin first derivative */
/* -------------------------------------------------------------------------- */
template <>
class Kelvin<3, 1> : protected Kelvin<3, 0> {
using parent = Kelvin<3, 0>;
protected:
static constexpr UInt dim = parent::dim;
public:
using parent::parent;
template <bool upper, typename ST>
inline Vector<Complex, dim>
applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <bool upper, typename ST>
inline Vector<Complex, dim>
applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
};
template <bool upper, typename ST>
inline Vector<Complex, Kelvin<3, 1>::dim>
Kelvin<3, 1>::applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp = f * computeD<upper>(q);
auto res = Kelvin<3, 0>::applyU0<false>(q, tmp);
res *= -sign(upper);
Vector<Real, dim> e3{{{0, 0, 1}}};
tmp.template mul<false>(f, e3);
res += Kelvin<3, 0>::applyU1<upper>(q, tmp);
res *= q.l2norm();
return res;
}
template <bool upper, typename ST>
inline Vector<Complex, Kelvin<3, 1>::dim>
Kelvin<3, 1>::applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp = f * computeD<upper>(q);
auto res = Kelvin<3, 0>::applyU1<upper>(q, tmp);
res *= -sign(upper) * q.l2norm();
return res;
}
/* -------------------------------------------------------------------------- */
/* Kelvin second derivative */
/* -------------------------------------------------------------------------- */
template <>
class Kelvin<3, 2> : protected Kelvin<3, 1> {
using parent = Kelvin<3, 1>;
static constexpr UInt dim = parent::dim;
public:
using parent::parent;
template <bool upper, typename ST>
inline Matrix<Complex, dim, dim>
applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <bool upper, typename ST>
inline Matrix<Complex, dim, dim>
applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
template <typename ST>
inline Matrix<Complex, dim, dim>
applyDiscontinuityTerm(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const;
};
template <bool upper, typename ST>
inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
Kelvin<3, 2>::applyU0(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp = Kelvin<3, 1>::applyU0<upper>(q, f);
tmp *= -sign(upper);
auto res = outer(computeD<upper>(q), tmp);
Vector<Real, dim> e3{{{0, 0, 1}}};
res += outer(e3, Kelvin<3, 1>::applyU1<upper>(q, f));
res *= -1. * q.l2norm(); // << -1 for grad_x
return res;
}
template <bool upper, typename ST>
inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
Kelvin<3, 2>::applyU1(VectorProxy<const Real, dim - 1>& q,
StaticMatrix<Complex, ST, dim, dim>& f) const {
auto tmp = Kelvin<3, 1>::applyU1<upper>(q, f);
auto res = outer(computeD<upper>(q), tmp);
res *= -1. * -sign(upper) * q.l2norm(); // << -1 for grad_x
return res;
}
template <typename ST>
inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
Kelvin<3, 2>::applyDiscontinuityTerm(
VectorProxy<const Real, dim - 1>& /*q*/,
StaticMatrix<Complex, ST, dim, dim>& f) const {
Vector<Real, dim> e3{{{0, 0, 1}}};
auto res = (f * e3);
res *= -this->b;
res(2) += 2. * f(2, 2);
e3(2) *= -1 / (this->mu * this->b); // << -1 for grad_x
return outer(e3, res);
}
/* -------------------------------------------------------------------------- */
/* Implementation */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <UInt order>
struct KelvinIntegrator {
static inline Real int_exp(Real a, Real z) { return std::exp(a * z) / a; }
static inline Real int_exp_z(Real a, Real z) {
const Real a_inv = 1. / a;
return std::exp(a * z) * (z - a_inv) * a_inv;
}
static inline Real int_exp_zz(Real a, Real z) {
const Real a_inv = 1. / a;
return std::exp(a * z) * (z * z - 2 * a_inv * z + 2 * a_inv * a_inv) *
a_inv;
}
template <bool upper>
static inline Real g0(Real z) {
return std::exp(-sign(upper) * z);
}
template <bool upper>
static inline Real int_g0(const thrust::pair<Real, Real>& bounds,
Real alpha) {
alpha *= -sign(upper);
return int_exp(alpha, bounds.second) - int_exp(alpha, bounds.first);
}
template <bool upper>
static inline Real int_g0_z(const thrust::pair<Real, Real>& bounds,
Real alpha) {
alpha *= -sign(upper);
return int_exp_z(alpha, bounds.second) - int_exp_z(alpha, bounds.first);
}
template <bool upper>
static inline Real g1(Real z) {
return g0<upper>(z) * z;
}
template <bool upper>
static inline Real int_g1(const thrust::pair<Real, Real>& bounds,
Real alpha) {
return alpha * int_g0_z<upper>(bounds, alpha);
}
template <bool upper>
static inline Real int_g1_z(const thrust::pair<Real, Real>& bounds,
Real alpha) {
return alpha * (int_exp_zz(-alpha * sign(upper), bounds.second) -
int_exp_zz(-alpha * sign(upper), bounds.first));
}
template <Int yj_xi, UInt shape_num, UInt dim_q, typename OutType,
typename InType, typename KelvinType>
static inline void integrate(OutType&& out, InType&& in, KelvinType&& kelvin,
VectorProxy<const Real, dim_q>& q_, Real dij,
Real dl);
};
// Integrate against a constant shape function
// template <>
// template <Int yj_xi, UInt shape_num, UInt dim_q, typename OutType,
// typename InType, typename KelvinType>
// inline void KelvinIntegrator<0>::integrate(OutType&& out, InType&& in,
// KelvinType&& kelvin,
// VectorProxy<const Real, dim_q>&
// q_,
// Real dij, Real dl) {
// const auto q = q_.l2norm();
// const auto alpha = q * dl, beta = q * dij;
// // Integrate on [-1, 0]
// const auto bounds_low = thrust::make_pair(-0.5, 0);
// constexpr bool domain_low = yj_xi > 0;
// auto g0_mul{kelvin.template applyU0<domain_low>(q_, in)};
// auto g1_mul{kelvin.template applyU1<domain_low>(q_, in)};
// decltype(g0_mul) res;
// res = g0_mul * (g0<domain_low>(beta) * int_g0<domain_low>(bounds_low,
// alpha));
// res +=
// g1_mul * (g1<domain_low>(beta) * int_g0<domain_low>(bounds_low, alpha)
// +
// g0<domain_low>(beta) * int_g1<domain_low>(bounds_low,
// alpha));
// // Integrate on [0, 1]
// const auto bounds_up = thrust::make_pair(0, 0.5);
// constexpr bool domain_up = yj_xi >= 0;
// // Avoid computing
// if (yj_xi == 0) {
// g0_mul = kelvin.template applyU0<true>(q_, in);
// g1_mul = kelvin.template applyU1<true>(q_, in);
// }
// res += g0_mul * (g0<domain_up>(beta) * int_g0<domain_up>(bounds_up,
// alpha));
// res += g1_mul * (g1<domain_up>(beta) * int_g0<domain_up>(bounds_up, alpha)
// +
// g0<domain_up>(beta) * int_g1<domain_up>(bounds_up,
// alpha));
// res *= dl;
// out += res;
// }
namespace detail {
template <UInt shape_num>
struct LinearShape {
using integrator = KelvinIntegrator<1>;
template <bool domain, typename T>
static inline T integrateG0(T& g0_mul, Real alpha, Real beta);
template <bool domain, typename T>
static inline T integrateG1(T& g1_mul, Real alpha, Real beta);
};
template <UInt shape_num>
template <bool domain, typename T>
inline T LinearShape<shape_num>::integrateG0(T& g0_mul, Real alpha,
Real beta) {
constexpr Int mul = (shape_num == 0) ? -1 : 1;
const auto bounds = thrust::make_pair(std::min(-mul, 0), std::max(-mul, 0));
return g0_mul * (integrator::g0<domain>(beta) *
(integrator::int_g0<domain>(bounds, alpha) +
mul * integrator::int_g0_z<domain>(bounds, alpha)));
}
template <UInt shape_num>
template <bool domain, typename T>
inline T LinearShape<shape_num>::integrateG1(T& g1_mul, Real alpha,
Real beta) {
constexpr Int mul = (shape_num == 0) ? -1 : 1;
const auto bounds = thrust::make_pair(std::min(-mul, 0), std::max(-mul, 0));
return g1_mul * (integrator::g1<domain>(beta) *
(integrator::int_g0<domain>(bounds, alpha) +
mul * integrator::int_g0_z<domain>(bounds, alpha)) +
integrator::g0<domain>(beta) *
(integrator::int_g1<domain>(bounds, alpha) +
mul * integrator::int_g1_z<domain>(bounds, alpha)));
}
}
// Integrate against a linear shape function
template <>
template <Int yj_xi, UInt shape_num, UInt dim_q, typename OutType,
typename InType, typename KelvinType>
inline void KelvinIntegrator<1>::integrate(OutType&& out, InType&& in,
KelvinType&& kelvin,
VectorProxy<const Real, dim_q>& q_,
Real dij, Real dl) {
const auto q = q_.l2norm();
const auto alpha = q * dl, beta = q * dij;
constexpr bool domain = (shape_num == 0) ? yj_xi >= 0 : yj_xi > 0;
auto g0_mul = kelvin.template applyU0<domain>(q_, in);
auto g1_mul = kelvin.template applyU1<domain>(q_, in);
decltype(g0_mul) res;
res = detail::LinearShape<shape_num>::template integrateG0<domain>(
g0_mul, alpha, beta);
res += detail::LinearShape<shape_num>::template integrateG1<domain>(
g1_mul, alpha, beta);
res *= dl;
out += res;
}
/* -------------------------------------------------------------------------- */
/// Class for the Boussinesq tensor
template <UInt dim, UInt derivative_order>
class Boussinesq;
template <>
class Boussinesq<3, 0> {
protected:
static constexpr UInt dim = 3;
public:
/// Constructor
Boussinesq(Real mu, Real nu)
: mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
template <typename ST>
inline Vector<Complex, dim>
applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
template <typename ST>
inline Vector<Complex, dim>
applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
protected:
const Real mu, nu, lambda;
};
template <typename ST>
inline Vector<Complex, Boussinesq<3, 0>::dim>
Boussinesq<3, 0>::applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// This guy has the property Bt(q) = B(-q). Since we know we recieve -q as an
// argument, we can just apply B without transpose
auto dplus = computeD<true>(q), dminus = computeD<false>(q), d2 = computeD2(q);
auto res = 2. * t;
res += ((1. - 2. * nu) * dminus.dot(t)) * dplus;
res += d2.dot(t) * d2;
res(2) -= t(2);
res *= (1. / (2. * mu * q.l2norm()));
return res;
}
template <typename ST>
inline Vector<Complex, Boussinesq<3, 0>::dim>
Boussinesq<3, 0>::applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// This guy does not have the property that Bt(q) = B(-q), so we have to apply
// Bt(q) for real. This makes us correct (in the line below) for the fact that
// we recieve -q as an argument.
auto dplus = computeD<false>(q);
auto res = (dplus.dot(t) / (2. * mu * q.l2norm())) * dplus;
return res;
}
/* -------------------------------------------------------------------------- */
template <>
class Boussinesq<3, 1> : protected Boussinesq<3, 0> {
protected:
using parent = Boussinesq<3, 0>;
static const UInt dim = parent::dim;
public:
using parent::parent;
template <typename ST>
inline Matrix<Complex, dim, dim>
applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
template <typename ST>
inline Matrix<Complex, dim, dim>
applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const;
};
template <typename ST>
inline Matrix<Complex, Boussinesq<3, 1>::dim, Boussinesq<3, 1>::dim>
Boussinesq<3, 1>::applyU0(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// taking into account K(-q) convolution
auto dbar = computeD<false>(q);
Vector<Real, dim> e3{{{0, 0, 1}}};
auto res = outer(dbar, Boussinesq<3, 0>::applyU0(t, q));
res *= -1.;
res += outer(e3, Boussinesq<3, 0>::applyU1(t, q));
res *= q.l2norm();
return res;
}
template <typename ST>
inline Matrix<Complex, Boussinesq<3, 1>::dim, Boussinesq<3, 1>::dim>
Boussinesq<3, 1>::applyU1(StaticVector<Complex, ST, dim>& t,
VectorProxy<const Real, dim - 1>& q) const {
// taking into account K(-q) convolution
auto dbar = computeD<false>(q);
auto res = outer(dbar, Boussinesq<3, 0>::applyU1(t, q));
res *= -q.l2norm();
return res;
}
} // namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif // INFLUENCE_HH
c4science · Help