3D EM based on Schur decomposition
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

EMSchur3DBase.h 14KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473
  1. // ===========================================================================
  2. //
  3. // Filename: EMSchur3DBase.h
  4. //
  5. // Created: 09/20/2013 04:35:57 PM
  6. // Compiler: Tested with g++, icpc, and MSVC 2010
  7. //
  8. // Author: Trevor Irons (ti)
  9. //
  10. // Organisation: University of Utah,
  11. // Colorado School of Mines
  12. // US Geological Survey
  13. //
  14. // Email: Trevor.Irons@utah.edu
  15. //
  16. // ===========================================================================
  17. /**
  18. @file
  19. @author Trevor Irons
  20. @date 09/20/2013
  21. @version $Id$
  22. **/
  23. #pragma once
  24. #ifndef EMSCHUR3DBASE_INC
  25. #define EMSCHUR3DBASE_INC
  26. #include <LemmaCore>
  27. #include "EMSchur3DConfig.h"
  28. #include <FDEM1D>
  29. //#include "LemmaObject.h"
  30. //#include "rectilineargrid.h"
  31. //#include "RectilinearGridVTKExporter.h"
  32. //#include "ASCIIParser.h"
  33. //#include "AEMSurvey.h"
  34. //#include "receiverpoints.h"
  35. //#include "layeredearthem.h"
  36. //#include "emearth1d.h"
  37. #include "timer.h"
  38. #include <Eigen/Sparse>
  39. //#define SPARSEMAJOR Eigen::ColMajor
  40. #define SPARSEMAJOR Eigen::RowMajor
  41. // SPARSEMAJOR;
  42. //typedef Eigen::RowMajor SPARSEMAJOR;
  43. //#include "bicgstab.h"
  44. // // Solvers
  45. // #ifdef HAVE_PASTIX
  46. // #include <Eigen/PaStiXSupport>
  47. // #endif
  48. //
  49. // #ifdef HAVE_METIS
  50. // #include <Eigen/MetisSupport>
  51. // #endif
  52. //
  53. // #ifdef HAVE_SUPERLU
  54. // #include <Eigen/SuperLUSupport>
  55. // #endif
  56. //
  57. // #ifdef HAVE_SUPERLUMT
  58. // #include <Eigen/SuperLUMTSupport>
  59. // #endif
  60. //
  61. // #ifdef HAVE_SPQR
  62. // #include <Eigen/SPQRSupport>
  63. // #endif
  64. // Cholmod Support won't compile typedef issue
  65. // #ifdef HAVE_CHOLMOD
  66. // #include <Eigen/CholmodSupport>
  67. // #endif
  68. //
  69. // // Cholmod Support won't compile
  70. // #ifdef HAVE_UMFPACK
  71. // #include <Eigen/UmfPackSupport>
  72. // #endif
  73. namespace Lemma {
  74. // Pardiso likes LOWER for Sym problem,
  75. // for non-symmetric, set both to 1
  76. #define UPPER 1 // 1=true, 0=false
  77. #define LOWER 1 // 1=true, 0=false
  78. enum SOLVER{ SPARSELU, SimplicialLLT, SimplicialLDLT, BiCGStab, SparseQR };
  79. /**
  80. \ingroup EMSchur3D
  81. \brief Provides a 3D solution to Maxwell's equations.
  82. \details 3D finite difference solution to maxwells equations
  83. using a SCHUR decomposition on a staggered grid.
  84. Performs a Schur decomposition on the vector scalar formulation of
  85. Maxwell's equations.
  86. \f[
  87. -\nabla^2 (\mathbf{A}) - \jmath \omega \mu \sigma \mathbf{A} - \mu \sigma \nabla (\phi) = - \mu \mathbf{J}_s
  88. \f]
  89. Which can be written in the functional form
  90. \f[ \begin{pmatrix}
  91. -\nabla^2 + \jmath \omega \mu \sigma & \mu \sigma \nabla \\
  92. \nabla \cdot & 0
  93. \end{pmatrix}
  94. \begin{pmatrix} \mathbf{A} \\ \phi \end{pmatrix}
  95. = \begin{pmatrix} \mathbf{s}_E \\ 0 \end{pmatrix}
  96. \f]
  97. Using the notation
  98. \f[ \begin{pmatrix}
  99. \mathbf{C} & \mathbf{B} \\
  100. \mathbf{D} & \mathbf{0}
  101. \end{pmatrix} \begin{pmatrix} \mathbf{A} \\ \phi \end{pmatrix} =
  102. \begin{pmatrix} \mathbf{s}_E \\ 0 \end{pmatrix}
  103. \f]
  104. Which is decomposed for seperate solutions to \f$ \mathbf{A}, \phi \f$ using a Schur decomposition
  105. \f[ \begin{matrix}
  106. \mathbf{D}\mathbf{C}^{-1}\mathbf{B} \phi & = \mathbf{D} \mathbf{C}^{-1} \mathbf{s}_E \\
  107. \mathbf{C}\mathbf{A} & = \mathbf{s}_E - \mathbf{G} \phi
  108. \end{matrix}
  109. \f]
  110. Where \f$ \mathbf{B} = \mu \sigma \mathbf{D}^T \f$. Additional algorithmic details may be found at
  111. @verbatim
  112. @inproceedings{doi:10.1190/segam2012-0896.1,
  113. author = {Trevor Irons and Yaoguo Li and Jason R. McKenna},
  114. title = {3D frequency-domain electromagnetics modeling using decoupled scalar and vector potentials},
  115. booktitle = {SEG Technical Program Expanded Abstracts 2012},
  116. chapter = {112},
  117. year = {2012},
  118. pages = {1-6},
  119. doi = {10.1190/segam2012-0896.1},
  120. URL = {http://library.seg.org/doi/abs/10.1190/segam2012-0896.1},
  121. eprint = {http://library.seg.org/doi/pdf/10.1190/segam2012-0896.1}
  122. }
  123. @endverbatim
  124. */
  125. //template< class Solver >
  126. class EMSchur3DBase : public LemmaObject {
  127. friend std::ostream &operator<<(std::ostream &stream,
  128. const EMSchur3DBase &ob);
  129. protected:
  130. //template<typename U>
  131. //friend class EMSchur3D;
  132. public:
  133. // ==================== LIFECYCLE =======================
  134. /** Default protected constructor, use New */
  135. explicit EMSchur3DBase ( const ctor_key& );
  136. /** Default protected constructor, use New */
  137. explicit EMSchur3DBase ( const YAML::Node& node, const ctor_key& );
  138. /** Default protected destructor, use Delete */
  139. virtual ~EMSchur3DBase ( );
  140. /**
  141. * Initialises antenna to contain no points, with no current
  142. * and no frequency. NumberOfTurns set to 1
  143. */
  144. static std::shared_ptr<EMSchur3DBase> NewSP();
  145. /*
  146. * Provides deep copy
  147. */
  148. //virtual std::shared_ptr<EMSchur3DBase> Clone() const ;
  149. /**
  150. * Uses YAML to serialize this object.
  151. * @return a YAML::Node
  152. */
  153. YAML::Node Serialize() const;
  154. /**
  155. * Constructs an object from a YAML::Node.
  156. */
  157. static std::shared_ptr<EMSchur3DBase> DeSerialize( const YAML::Node& node );
  158. // ==================== OPERATORS =======================
  159. // ==================== OPERATIONS =======================
  160. /* Performs the solution
  161. * This is thread safe. TODO, but where should the results go? Just to file?
  162. * Who does the parsing? Actually I think this method is the wrong place to talk
  163. * about that. This is just a big red button. The details should be worked out in private
  164. * methods, that this could well call. Still though, where should the damn results go. What if
  165. * someone wants to use them *right now*, and not go through file IO. This is a good cause for
  166. * fixing the model class. So the result will be a final RectGrid (or Model) where the results live.
  167. * THEN users can call a seperate WriteToVTK or whatever based on that.
  168. */
  169. void Solve( );
  170. // ==================== ACCESS =======================
  171. /** Sets the RectilinearGrid to use
  172. * @param[in] Grid is a pointer to the Grid to be used.
  173. */
  174. void SetRectilinearGrid( std::shared_ptr<RectilinearGrid> Grid);
  175. /** Sets the RectilinearGrid to use
  176. * @param[in] Grid is a pointer to the Grid to be used.
  177. */
  178. void SetRectilinearGrid( vtkRectilinearGrid* vtkGrid );
  179. /** Sets the RectilinearGrid to use
  180. * @param[in] Grid is a pointer to the Grid to be used.
  181. */
  182. void SetAEMSurvey( std::shared_ptr<AEMSurvey> Survey);
  183. /** Sets the prefix for results files (.log and .vtr) the source fiducial is added as well
  184. */
  185. void SetResFileName(const std::string& filename);
  186. /** Sets the solver to use to invert the C matrix. This is done many times. If you are reusing the same matrix
  187. for numerous simulations, it may be benefitial to use the direct (SPARSELU) solver. For one off calculations the BiCGSTAB
  188. is a good choice. Default is SPARSELU.
  189. */
  190. void SetCSolver(const SOLVER& CSOLVER);
  191. /**
  192. * Sets the LayredEarthEM model that will be used for the primary field calculation as well as deterimining the
  193. * bacground conductivity file.
  194. @ @param[in] LayModel is a pointer to the LayeredEarthEM model to use.
  195. */
  196. void SetLayeredEarthEM( std::shared_ptr<LayeredEarthEM> LayModel );
  197. /**
  198. * Loads a MeshTools conductivity model.
  199. * @param[in] fname is the file to load.
  200. */
  201. void LoadMeshToolsConductivityModel( const std::string& fname );
  202. /**
  203. * Writes out results to into a vtkRectilinearGrid file
  204. * @param[in] fname is the file name that is created, the .vtr suffix is added.
  205. */
  206. void WriteVTKResults( const std::string& fname, Eigen::Ref<VectorXcr> A,
  207. Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0, Eigen::Ref<VectorXcr> E,
  208. Eigen::Ref<VectorXcr> Phi, Eigen::Ref<VectorXcr> ADiv, Eigen::Ref<VectorXcr> ADiv2,
  209. Eigen::Ref<VectorXcr> B
  210. );
  211. // ==================== INQUIRY =======================
  212. /** Returns the name of the underlying class, similiar to Python's type */
  213. virtual std::string GetName() const {
  214. return this->CName;
  215. }
  216. protected:
  217. // ==================== LIFECYCLE =======================
  218. //private:
  219. /** Builds the C matrix */
  220. void BuildC( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw);
  221. /* Builds the C matrix as a block real system. Benefits of this are the availability of an
  222. * LDL^T decomposition. Also, as complex number in C++ are templates and will necessarily have
  223. * real and imaginary parts, this formulation will have a reduced cost, due to less computations
  224. * with the zero valued imaginary parts (off-diagonals)
  225. * The \f$C \in I^3\f$ matrix is instead written as
  226. * [ C_r C_i ] [ x_r ] = [ b_r ]
  227. * [ C_i -C_r ] [ x_i ] [ b_i ]
  228. */
  229. //void BuildCReal( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw);
  230. /** Builds the C matrix */
  231. void BuildCPreconditioner( const int& iw);
  232. /** Builds the C matrix */
  233. virtual void BuildCDirectSolver( )=0;
  234. /** Fills the actual points on the grid that 1D primary field calculations need to be made.
  235. @todo This is a little stupid as all threads share the same points. Stupid in that right now
  236. this is done for every calculation. Not a huge amount of time is spent here, I suppose some extra memory
  237. though. We need to extend
  238. EmEARTH to be able to input a grid so that points are not explicitly needed like
  239. this. This requires some care as calcs are made on faces.
  240. Alternatively, the bins function of FieldPoints could be extended quite easily.
  241. This may be the way to do this.
  242. */
  243. void FillPoints( std::shared_ptr<FieldPoints> Points );
  244. /** Builds D/G
  245. */
  246. void BuildD( );
  247. /** Used to manage tradititional C 3D array */
  248. template <typename T>
  249. void Allocate3DScalar(T ***&Array, T init) {
  250. Array = new T**[nx];
  251. for (int ix=0; ix<nx; ix++){
  252. Array[ix] = new T*[ny];
  253. for (int iy=0; iy<ny; iy++){
  254. Array[ix][iy] = new T[nz];
  255. for (int iz=0; iz<nz; iz++) Array[ix][iy][iz] = init;
  256. }
  257. }
  258. return;
  259. }
  260. /** Used to manage tradititional C 3D array */
  261. template <typename T>
  262. void Delete3DScalar(T ***&Array) {
  263. for (int ix=0; ix<nx; ix++){
  264. for (int iy=0; iy<ny; iy++){
  265. delete [] Array[ix][iy];
  266. }
  267. delete [] Array[ix];
  268. }
  269. delete [] Array;
  270. Array = NULL;
  271. return;
  272. }
  273. /**
  274. * This is called just before solve and gets all shared objects ready to go
  275. */
  276. void Setup( );
  277. /** Solves a single source problem. This method is thread safe.
  278. * @param[in] Source is the source term for generating primary fields
  279. * @param[in] isource is the source index
  280. */
  281. virtual void SolveSource( std::shared_ptr< DipoleSource > Source , const int& isource)=0;
  282. /** Computes the primary field
  283. */
  284. void PrimaryField( std::shared_ptr<DipoleSource> Source, std::shared_ptr<FieldPoints> dpoint);
  285. /**
  286. * Fills the vectors that are called in
  287. */
  288. void FillSourceTerms( Eigen::Ref<VectorXcr> ms,
  289. Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0,
  290. std::shared_ptr<FieldPoints> dpoint, const Real& omega );
  291. /** Computes the curl of A on the staggered grid
  292. */
  293. VectorXcr StaggeredGridCurl(Eigen::Ref<VectorXcr> A);
  294. // ==================== DATA MEMBERS =========================
  295. /** Grid over which operators are active */
  296. std::shared_ptr<RectilinearGrid> Grid;
  297. /* Used to help dump results */
  298. //std::shared_ptr<RectilinearGridVTKExporter> VTKGridExporter;
  299. /** Class containing information about the AEM survey */
  300. std::shared_ptr<AEMSurvey> Survey;
  301. std::shared_ptr<LayeredEarthEM> LayModel;
  302. /** Matrix objects representing discrete operators
  303. */
  304. Eigen::SparseMatrix<Complex, SPARSEMAJOR>* Cvec;
  305. Eigen::SparseMatrix<Complex, SPARSEMAJOR> D;
  306. /** Squeezed matrices for reduced phi grid
  307. */
  308. Eigen::SparseMatrix<Complex, SPARSEMAJOR>* Cvec_s;
  309. Eigen::SparseMatrix<Complex, SPARSEMAJOR> D_s;
  310. /** number of cells in x, set when RectilinearGrid is attached */
  311. int nx;
  312. /** number of cells in y, set when RectilinearGrid is attached */
  313. int ny;
  314. /** number of cells in z, set when RectilinearGrid is attached */
  315. int nz;
  316. /** number of fields/faces in x, unwrapped x */
  317. int unx;
  318. /** number of fields/faces in y, unwrapped y */
  319. int uny;
  320. /** number of fields/faces in z, unwrapped z */
  321. int unz;
  322. /** number of cell centres, unwrapped scalars */
  323. int uns;
  324. /** name for log files and VTK files */
  325. std::string ResFile;
  326. /** frequency of source */
  327. VectorXr Omegas;
  328. /** Conductivity model */
  329. //Complex ***sigma_jwe;
  330. Real ***sigma;
  331. /** Conductivity model minus reference model */
  332. //Complex ***sigmap;
  333. Real ***sigmap;
  334. /** rectilinear grid spacing in x */
  335. VectorXr hx;
  336. /** rectilinear grid spacing in y */
  337. VectorXr hy;
  338. /** rectilinear grid spacing in z */
  339. VectorXr hz;
  340. /** inverse of hx */
  341. VectorXr ihx;
  342. /** inverse of hx squared */
  343. VectorXr ihx2;
  344. /** inverse of hy */
  345. VectorXr ihy;
  346. /** inverse of hy squared */
  347. VectorXr ihy2;
  348. /** inverse of hz */
  349. VectorXr ihz;
  350. /** inverse of hz squared */
  351. VectorXr ihz2;
  352. /** Marker for air cells */
  353. VectorXi MAC;
  354. /** Marker for air cells */
  355. std::vector<int> idx;
  356. /** Dirichlet on low X */
  357. bool DirichletXLO = true;
  358. //bool DirichletXLO = false;
  359. bool DirichletXHI = true;
  360. //bool DirichletXHI = false;
  361. bool DirichletYLO = true;
  362. //bool DirichletYLO = false;
  363. bool DirichletYHI = true;
  364. //bool DirichletYHI = false;
  365. bool DirichletZLO = true;
  366. //bool DirichletZLO = false;
  367. bool DirichletZHI = true;
  368. //bool DirichletZHI = false;
  369. private:
  370. static constexpr auto CName = "EMSchur3DBase";
  371. }; // ----- end of class EMSchur3DBase -----
  372. }
  373. #endif // ----- #ifndef EMSCHUR3BASE_INC -----