Main Lemma Repository
Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.

KernelEM1DReflSpec.h 9.9KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author Trevor Irons
  8. @date 05/18/2012
  9. **/
  10. #ifndef KERNELEM1DREFLSPEC_INC
  11. #define KERNELEM1DREFLSPEC_INC
  12. #include "DipoleSource.h"
  13. #include "KernelEM1DReflBase.h"
  14. #include "LayeredEarthEM.h"
  15. namespace Lemma {
  16. // forward declaration for friend
  17. //template<EMMODE Mode, int Ikernel, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  18. //class KernelEM1DSpec;
  19. // ===================================================================
  20. // Class: KernelEM1DReflSpec
  21. /**
  22. @class
  23. \brief Specialized version of KernelEM1DReflBase
  24. \details Through use of template specialisations, this KernelEm1D
  25. class delivers much better performance. This class is internal
  26. to Lemma, you should never need to instantiate it. The constructors
  27. are public to allow make_shared. Additonally, this class is not
  28. serializable.
  29. */
  30. // ===================================================================
  31. template<EMMODE Mode, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  32. class KernelEM1DReflSpec : public KernelEM1DReflBase {
  33. // TODO can we use the manager's key instead to lock for that?
  34. struct ctor_key{};
  35. public:
  36. //template<EMMODE Mode2, int Ikernel2, DIPOLE_LOCATION Isource2, DIPOLE_LOCATION Irecv2>
  37. //friend class KernelEM1DSpec;
  38. friend class KernelEM1DManager;
  39. // ==================== LIFECYCLE =======================
  40. /// Default locked constructor.
  41. explicit KernelEM1DReflSpec ( const ctor_key& ) : KernelEM1DReflBase( ) {
  42. }
  43. /// Default protected constructor.
  44. ~KernelEM1DReflSpec () {
  45. }
  46. static std::shared_ptr<KernelEM1DReflSpec<Mode, Isource, Irecv> > NewSP() {
  47. return std::make_shared<KernelEM1DReflSpec<Mode, Isource, Irecv> >( ctor_key() );
  48. }
  49. // ==================== OPERATORS =======================
  50. // ==================== OPERATIONS =======================
  51. // ==================== ACCESS =======================
  52. // ==================== INQUIRY =======================
  53. virtual std::string GetName() {
  54. return CName;
  55. }
  56. protected:
  57. private:
  58. // ==================== LIFECYCLE =======================
  59. // ==================== OPERATIONS =======================
  60. /** Computes reflection coefficients. Depending on the
  61. * specialisation, this will either be TM or TE mode.
  62. */
  63. void ComputeReflectionCoeffs(const Real &lambda);
  64. /* Computes reflection coefficients. Depending on the
  65. * specialisation, this will either be TM or TE mode. This method
  66. * stores previous results in a struct. For a given index, and
  67. * lambda, the result will be the same. Turned out to be of limited utility.
  68. */
  69. //void ComputeReflectionCoeffs(const Real &lambda, const int& idx);
  70. /** Precomputes expensive calculations that are reused by insances
  71. * of KernelEM1DSpec in the calculation of Related potentials. This
  72. * method is specialised based on template parameters
  73. */
  74. void PreComputePotentialTerms();
  75. /*
  76. * Sets the cache in CACHE to use. Somewhat expensive, avoid calling in tight loops
  77. */
  78. //void SetTCache(const Real& rho0);
  79. // ==================== DATA MEMBERS =========================
  80. static constexpr auto CName = "KernelEM1DSpec";
  81. }; // ----- end of class KernelEM1DReflSpec -----
  82. //template<EMMODE Mode, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  83. //std::unordered_map<Real, cache> KernelEM1DReflSpec<Mode, Isource, Irecv>::CACHE;
  84. //template<EMMODE Mode, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  85. //cache* KernelEM1DReflSpec<Mode, Isource, Irecv>::tcache;
  86. ///////////////////////////////////////////////
  87. // Declarations of specialisations
  88. template<>
  89. void KernelEM1DReflSpec<TM, INAIR, INAIR>::ComputeReflectionCoeffs(const Real& lambda);
  90. template<>
  91. void KernelEM1DReflSpec<TE, INAIR, INAIR>::ComputeReflectionCoeffs(const Real& lambda);
  92. template<>
  93. void KernelEM1DReflSpec<TM, INAIR, INGROUND>::ComputeReflectionCoeffs(const Real& lambda);
  94. template<>
  95. void KernelEM1DReflSpec<TE, INAIR, INGROUND>::ComputeReflectionCoeffs(const Real& lambda);
  96. template<>
  97. void KernelEM1DReflSpec<TM, INAIR, INAIR>::PreComputePotentialTerms( );
  98. template<>
  99. void KernelEM1DReflSpec<TE, INAIR, INAIR>::PreComputePotentialTerms( );
  100. template<>
  101. void KernelEM1DReflSpec<TM, INAIR, INGROUND>::PreComputePotentialTerms( );
  102. template<>
  103. void KernelEM1DReflSpec<TE, INAIR, INGROUND>::PreComputePotentialTerms( );
  104. ///////////////////////////////////////////////
  105. // Default mode definitions
  106. template<EMMODE Mode, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  107. void KernelEM1DReflSpec<Mode, Isource, Irecv>::ComputeReflectionCoeffs(const Real& lambda) {
  108. static bool called = false;
  109. if (!called) {
  110. std::cout << "unspecialised Reflection function KernelEM1DReflSpec<"
  111. << Mode << ", " << Isource << ", "
  112. << Irecv << " >::ComputeReflectionCoeffs( const Real& lambda ) --> SLOW PERFORMANCE EXPECTED\n";
  113. called = true;
  114. }
  115. rams = lambda*lambda;
  116. //////////////////////////////////////////
  117. // Compute TEM stuff
  118. // This call to sqrt takes ~ 15% of execution time
  119. u = (rams-kk.array()).sqrt();
  120. uk = u(lays);
  121. um = u(layr);
  122. switch (Mode) {
  123. // TM mode
  124. case (TM):
  125. Zyu(1) = -u(0)/yh(0);
  126. Zyi = u.array() / yh.array();
  127. break;
  128. // TE mode
  129. case (TE):
  130. Zyu(1) = -u(0)/zh(0);
  131. Zyi = u.array() / zh.array();
  132. break;
  133. default:
  134. throw 11; //IllegalMode(this);
  135. }
  136. Zyd.tail<1>() = Zyi.tail<1>();
  137. for (int ilay=1; ilay<nlay-1; ++ilay) {
  138. cf(ilay) =
  139. std::exp(-(Real)(2.)*u(ilay)*LayerThickness(ilay));
  140. th(ilay) = ((Real)(1.)-cf(ilay)) / ((Real)(1.)+cf(ilay));
  141. }
  142. // Can't use blocks, b/c recursive
  143. for (int N=1; N<lays; ++N){
  144. Zyu(N+1)=Zyi(N)*(Zyu(N)-Zyi(N)*th(N)) /
  145. (Zyi(N)-Zyu(N)*th(N)) ;
  146. }
  147. int ne = nlay-2;
  148. for (int N=ne; N >=lays+1; --N) {
  149. Zyd(N) = Zyi(N)*(Zyd(N+1)+Zyi(N)*th(N)) /
  150. (Zyi(N)+Zyd(N+1)*th(N)) ;
  151. }
  152. rtd(nlay-1) = 0;
  153. if (layr < lays) {
  154. // Receiver above source layer
  155. int ls = layr;
  156. if (ls == 0) {
  157. ls = layr+1;
  158. }
  159. for (int N=ls; N<=lays; ++N) {
  160. rtu(N)= (Zyi(N)+Zyu(N)) /
  161. (Zyi(N)-Zyu(N)) ;
  162. }
  163. if (lays < nlay-1) {
  164. rtd(lays) = (Zyi(lays)-Zyd(lays+1)) /
  165. (Zyi(lays)+Zyd(lays+1)) ;
  166. }
  167. } else {
  168. // RECEIVER IN OR BELOW THE SOURCE LAYER
  169. if (lays == layr) { // Rx In source Layer
  170. if (layr == 0) {
  171. rtd(0) = (Zyu(1)+Zyd(1)) /
  172. (Zyu(1)-Zyd(1)) ;
  173. } else if (layr == nlay-1) {
  174. rtu(nlay-1) = (Zyi(nlay-1)+Zyu(nlay-1)) /
  175. (Zyi(nlay-1)-Zyu(nlay-1)) ;
  176. } else {
  177. rtu(layr) = (Zyi(layr)+Zyu(layr)) /
  178. (Zyi(layr)-Zyu(layr)) ;
  179. rtd(layr) = (Zyi(layr)-Zyd(layr+1)) /
  180. (Zyi(layr)+Zyd(layr+1)) ;
  181. }
  182. } else { // receiver below source layer
  183. if (lays == 0) {
  184. rtd(0) = (Zyu(1)+Zyd(1)) /
  185. (Zyu(1)-Zyd(1)) ;
  186. } else {
  187. rtu(lays) = (Zyi(lays)+Zyu(lays)) /
  188. (Zyi(lays)-Zyu(lays)) ;
  189. }
  190. int le = layr;
  191. if (le == nlay-1) --le;
  192. int ls = lays;
  193. if (lays == 0 ) ++ls;
  194. // TODO use blocks to vectorize maybe?
  195. // This works but gives same to slightly worse
  196. // performance as loop.
  197. // int nn = le-ls+1;
  198. // rtd.segment(ls, nn) =
  199. // (Zyi.segment(ls , nn).array() -
  200. // Zyd.segment(ls+1, nn).array()).array() /
  201. // (Zyi.segment(ls , nn).array() +
  202. // Zyd.segment(ls+1, nn).array()).array() ;
  203. for (int N=ls; N<=le; ++N) {
  204. rtd(N) = (Zyi(N)-Zyd(N+1)) /
  205. (Zyi(N)+Zyd(N+1)) ;
  206. }
  207. }
  208. } // End in or below source layer
  209. return;
  210. }
  211. template<EMMODE Mode, DIPOLE_LOCATION Isource, DIPOLE_LOCATION Irecv>
  212. void KernelEM1DReflSpec<Mode, Isource, Irecv>::PreComputePotentialTerms( ) {
  213. static bool called = false;
  214. if (!called) {
  215. std::cerr << "unspecialised function KernelEM1DReflSpec<"
  216. << Mode << ", " << Isource << ", "
  217. << Irecv << " >::PreComputePotentialTerms\n";
  218. called = true;
  219. }
  220. }
  221. } // ----- end of Lemma name -----
  222. #endif // ----- #ifndef KERNELEM1DREFLSPEC_INC -----