3D EM based on Schur decomposition
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

EMSchur3D.cpp 3.9KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 01/04/2018 04:34:54 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email tirons@egi.utah.edu
  14. * @copyright Copyright (c) 2018, University of Utah
  15. * @copyright Copyright (c) 2018, Lemma Software, LLC
  16. */
  17. #include <LemmaCore>
  18. #include "RectilinearGridReader.h"
  19. //#include "AEMSurveyReader.h"
  20. #include "EMSchur3D.h"
  21. //#include "LayeredEarthEMReader.h"
  22. //#include "CSymSimplicialCholesky.h"
  23. using namespace Lemma;
  24. int main( int argc, char** argv ) {
  25. // BiCGSTAB Diagonal preconditioner
  26. auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
  27. if (argc < 3) {
  28. std::cout << "EMSchur3D <rgrid> <1dmod> <aemsurvey> " << std::endl;
  29. exit( EXIT_SUCCESS );
  30. }
  31. ///////////////////////////////////////////////////
  32. // Read in Grid
  33. auto GridRead = RectilinearGridReader::NewSP();
  34. try {
  35. GridRead->ReadASCIIGridFile( argv[1] );
  36. } catch (std::exception& e) {
  37. std::cout << "Caught an error opening ASCII Grid file: ";
  38. std::cout << e.what() << std::endl;
  39. std::cout << "enter filename or 0 to exit programme\n";
  40. std::string inp;
  41. std::cin >> inp;
  42. if (inp != "0")
  43. GridRead->ReadASCIIGridFile( inp.c_str() );
  44. else exit(EXIT_FAILURE);
  45. }
  46. //////////////////////////////////////////////////
  47. // Read in Layered earth, for backgound model
  48. <<<<<<< HEAD
  49. // auto LayEarthRead = LayeredEarthEMReader::NewSP();
  50. // try {
  51. // LayEarthRead->ReadASCIIInputFile( argv[2] );
  52. // } catch (std::exception& e) {
  53. // std::cout << "Caught an error opening ASCII Layered Earth file: ";
  54. // std::cout << e.what() << std::endl;
  55. // std::cout << "enter filename or 0 to exit programme\n";
  56. // std::string inp;
  57. // std::cin >> inp;
  58. // if (inp != "0")
  59. // GridRead->ReadASCIIGridFile( inp.c_str() );
  60. // else exit(EXIT_FAILURE);
  61. // }
  62. =======
  63. /*
  64. auto LayEarthRead = LayeredEarthEMReader::NewSP();
  65. try {
  66. LayEarthRead->ReadASCIIInputFile( argv[2] );
  67. } catch (std::exception& e) {
  68. std::cout << "Caught an error opening ASCII Layered Earth file: ";
  69. std::cout << e.what() << std::endl;
  70. std::cout << "enter filename or 0 to exit programme\n";
  71. std::string inp;
  72. std::cin >> inp;
  73. if (inp != "0")
  74. GridRead->ReadASCIIGridFile( inp.c_str() );
  75. else exit(EXIT_FAILURE);
  76. }
  77. */
  78. >>>>>>> bd2e3fa6b0dbb11555df9d4e7a7fe68e29b61ffd
  79. /*
  80. //////////////////////////////////////////////////
  81. // Read in source specification
  82. auto AEMRead = AEMSurveyReader::NewSP();
  83. try{
  84. AEMRead->ReadASCIIAEMFile(argv[4]);
  85. } catch (std::exception& e) {
  86. std::cout << "Caught an error opening ASCII AEM file: ";
  87. std::cout << e.what() << std::endl;
  88. std::cout << "enter filename or 0 to exit programme\n";
  89. std::string inp;
  90. std::cin >> inp;
  91. if (inp != "0")
  92. AEMRead->ReadASCIIAEMFile( inp.c_str() );
  93. else exit(EXIT_FAILURE);
  94. }
  95. */
  96. /*
  97. //////////////////////////////////////////////////
  98. // And solve
  99. EM3D->SetRectilinearGrid(GridRead->GetGrid());
  100. EM3D->SetLayeredEarthEM( LayEarthRead->GetLayeredEarth() );
  101. EM3D->SetAEMSurvey( AEMRead->GetSurvey() );
  102. EM3D->LoadMeshToolsConductivityModel( argv[3] );
  103. EM3D->SetResFileName(argv[5]);
  104. //EM3D->SetCSolver( Lemma::SPARSELU ); // Lemma::BiCGSTAB );
  105. EM3D->Solve();
  106. */
  107. std::cout << *EM3D;
  108. exit(EXIT_SUCCESS);
  109. }