Lemma is an Electromagnetics API
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

Hantenna.cpp 8.8KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262
  1. // ===========================================================================
  2. //
  3. // Filename: hantenna.cpp
  4. //
  5. // Created: 10/07/2010 08:57:04 AM
  6. // Modified: 11 April 2018
  7. // Compiler: Tested with g++, icpc, and MSVC 2017
  8. //
  9. // Author: Trevor Irons (ti)
  10. //
  11. // Copyright (C) 2012,2018 Trevor Irons
  12. //
  13. // Organisation: Lemma Software
  14. //
  15. // Email: Trevor.Irons@lemmasoftware.org
  16. //
  17. // ===========================================================================
  18. /**
  19. @file
  20. @author Trevor Irons
  21. @date 10/07/2010
  22. $Id$
  23. **/
  24. #include "LemmaCore"
  25. #include "FDEM1D"
  26. #include "timer.h"
  27. #if defined(__clang__)
  28. /* Clang/LLVM. ---------------------------------------------- */
  29. const char* compiler = "clang";
  30. const char* ver = __VERSION__;
  31. #elif defined(__ICC) || defined(__INTEL_COMPILER)
  32. /* Intel ICC/ICPC. ------------------------------------------ */
  33. const char* compiler = "icpc";
  34. const char* ver = __VERSION__;
  35. #elif defined(__GNUC__) || defined(__GNUG__)
  36. /* GNU GCC/G++. --------------------------------------------- */
  37. const char* compiler = "gcc (GCC) ";// __VERSION__;
  38. const char* ver = __VERSION__;
  39. #elif defined(_MSC_VER)
  40. /* Microsoft Visual Studio. --------------------------------- */
  41. const char* compiler = "msvc ";
  42. const int ver = _MSC_FULL_VER;
  43. #elif defined(__PGI)
  44. /* Portland Group PGCC/PGCPP. ------------------------------- */
  45. const char* compiler = "pgc";
  46. #endif
  47. using namespace Lemma;
  48. std::vector<Real> readinpfile(const std::string& fname);
  49. std::vector<std::string> readinpfile2(const std::string& fname);
  50. int main(int argc, char** argv) {
  51. const char *buildString = __DATE__ ", " __TIME__;
  52. std::cout
  53. << "===========================================================================\n"
  54. << "Lemma " << LEMMA_VERSION << "\n"
  55. << "[" << compiler << " " << ver << " " << buildString << "]\n"
  56. << "This program is part of Lemma, a geophysical modelling and inversion API. \n"
  57. << " This Source Code Form is subject to the terms of the Mozilla Public\n"
  58. << " License, v. 2.0. If a copy of the MPL was not distributed with this\n"
  59. << " file, You can obtain one at http://mozilla.org/MPL/2.0/. \n"
  60. << "Copyright (C) 2018 Lemma Software \n"
  61. << "More information may be found at: https://lemmasoftware.org\n"
  62. << " info@lemmasoftware.org\n"
  63. << "===========================================================================\n\n"
  64. << "Hantenna calculates the harmonic H field from polygonal wire loop sources\n";
  65. if (argc < 5) {
  66. std::cout << "usage: hantenna.exe trans.inp cond.inp points.inp config.inp \n";
  67. exit(0);
  68. }
  69. #ifdef LEMMAUSEOMP
  70. std::cout << "OpenMP is using " << omp_get_max_threads() << " threads" << std::endl;
  71. #endif
  72. std::vector<Real> Trans = readinpfile(std::string(argv[1]));
  73. std::vector<Real> CondMod = readinpfile(std::string(argv[2]));
  74. std::vector<Real> Points = readinpfile(std::string(argv[3]));
  75. std::vector<std::string> config = readinpfile2(std::string(argv[4]));
  76. //////////////////////////////////////
  77. // Define transmitter
  78. auto trans = PolygonalWireAntenna::NewSP();
  79. trans->SetNumberOfPoints((int)(Trans[0]));
  80. int ip=1;
  81. for ( ; ip<=(int)(Trans[0])*2; ip+=2) {
  82. trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], -1e-3));
  83. //trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], 50.));
  84. }
  85. trans->SetNumberOfFrequencies(1);
  86. trans->SetFrequency(0, Trans[ip]);
  87. trans->SetCurrent(Trans[ip+1]);
  88. trans->SetMinDipoleRatio(atof(config[1].c_str()));
  89. trans->SetMinDipoleMoment(atof(config[2].c_str()));
  90. trans->SetMaxDipoleMoment(atof(config[3].c_str()));
  91. /////////////////////////////////////
  92. // Field calculations
  93. auto receivers = FieldPoints::NewSP();
  94. int nx = (int)Points[0];
  95. int ny = (int)Points[1];
  96. int nz = (int)Points[2];
  97. Real ox = Points[3];
  98. Real oy = Points[4];
  99. Real oz = Points[5];
  100. Vector3r loc;
  101. VectorXr dx(nx-1);
  102. VectorXr dy(ny-1);
  103. VectorXr dz(nz-1);
  104. ip = 6;
  105. int ir = 0;
  106. for ( ; ip <6+nx-1; ++ip) {
  107. dx[ir] = Points[ip];
  108. ++ir;
  109. }
  110. ir = 0;
  111. for ( ; ip <6+ny-1+nx-1; ++ip) {
  112. dy[ir] = Points[ip];
  113. ++ir;
  114. }
  115. ir = 0;
  116. for ( ; ip <6+nz-1+ny-1+nx-1; ++ip) {
  117. dz[ir] = Points[ip];
  118. ++ir;
  119. }
  120. receivers->SetNumberOfPoints(nx*ny*nz);
  121. ir = 0;
  122. Real pz = oz;
  123. for (int iz=0; iz<nz; ++iz) {
  124. Real py = oy;
  125. for (int iy=0; iy<ny; ++iy) {
  126. Real px = ox;
  127. for (int ix=0; ix<nx; ++ix) {
  128. loc << px, py, pz;
  129. receivers->SetLocation(ir, loc);
  130. if (ix < nx-1) px += dx[ix];
  131. ++ ir;
  132. }
  133. if (iy<ny-1) py += dy[iy];
  134. }
  135. if (iz<nz-1) pz += dz[iz];
  136. }
  137. ////////////////////////////////////
  138. // Define model
  139. auto earth = LayeredEarthEM::NewSP();
  140. VectorXcr sigma;
  141. VectorXr thick;
  142. earth->SetNumberOfLayers(static_cast<int>(CondMod[0])+1);
  143. sigma.resize(static_cast<int>(CondMod[0])+1); sigma(0) = 0; // airlayer
  144. thick.resize(static_cast<int>(CondMod[0])-1);
  145. int ilay=1;
  146. for ( ; ilay/2<CondMod[0]-1; ilay+=2) {
  147. sigma(ilay/2+1) = 1./CondMod[ilay];
  148. thick(ilay/2) = CondMod[ilay+1];
  149. }
  150. sigma(ilay/2+1) = 1./ CondMod[ilay];
  151. earth->SetLayerConductivity(sigma);
  152. if (thick.size() > 0) earth->SetLayerThickness(thick);
  153. auto EmEarth = EMEarth1D::NewSP();
  154. EmEarth->AttachWireAntenna(trans);
  155. EmEarth->AttachLayeredEarthEM(earth);
  156. EmEarth->AttachFieldPoints(receivers);
  157. EmEarth->SetFieldsToCalculate(H);
  158. EmEarth->SetHankelTransformMethod(string2Enum<HANKELTRANSFORMTYPE>(config[0]));
  159. ///////////////////////////////////////////////
  160. // Keep track of time
  161. jsw_timer timer;
  162. timer.begin();
  163. clock_t launch = clock();
  164. EmEarth->CalculateWireAntennaFields(true); // true=status bar
  165. Real paTime = timer.end();
  166. std::cout << "\n\n===========================================\ncalc. real time: " << paTime/60. << "\t[m]\n";
  167. std::cout << "calc. user time: " << (clock()-launch)/CLOCKS_PER_SEC/60. << "\t[CPU m]"
  168. << std::endl;
  169. ////////////////////////////////////
  170. // Report
  171. std::fstream hrep("hfield.yaml", std::ios::out);
  172. std::fstream hreal("hfield.dat", std::ios::out);
  173. hrep << *EmEarth << std::endl;
  174. hrep.close();
  175. //hreal << *trans << std::endl;
  176. //hreal << *earth << std::endl;
  177. hreal << "// Right hand coordinate system, z is positive down\n";
  178. hreal << "// x[m]\ty[m]\tz[m]\tRe(Hx[A/m])\tRe(Hy[A/m])\tRe(Hz[A/m])\tIm(Hx)\tIm(Hy)\tIm(Hz)\n";
  179. hreal.precision(8);
  180. int i=0;
  181. for (int iz=0; iz<nz; ++iz) {
  182. for (int iy=0; iy<ny; ++iy) {
  183. for (int ix=0; ix<nx; ++ix) {
  184. hreal << receivers->GetLocation(i).transpose() << "\t";
  185. //hreal << receivers->GetHfield(0, i).transpose() << "\n"; // ( complex, notation )
  186. hreal << receivers->GetHfield(0, i).transpose().real() << "\t";
  187. hreal << receivers->GetHfield(0, i).transpose().imag() << "\n";
  188. ++i;
  189. }
  190. }
  191. }
  192. hreal.close();
  193. // Clean up
  194. // report timings
  195. #ifdef LEMMAUSEOMP
  196. std::ofstream outfile;
  197. outfile.open("timings.csv", std::ios_base::app);
  198. outfile << compiler << "," << ver << "," << buildString << "," << config[0] << "," << omp_get_max_threads() << "," << paTime/60. << "\n";
  199. #endif
  200. }
  201. std::vector<Real> readinpfile(const std::string& fname) {
  202. std::string buf;
  203. char dump[255];
  204. std::vector<Real> vals;
  205. std::fstream input(fname.c_str(), std::ios::in);
  206. if (input.fail()) {
  207. std::cerr << "Input file " << fname << " failed to open\n";
  208. exit(EXIT_FAILURE);
  209. }
  210. while (input >> buf) {
  211. if (buf.substr(0,2) == "//") {
  212. input.getline(dump, 255);
  213. } else {
  214. vals.push_back( atof(buf.c_str() ));
  215. }
  216. }
  217. return vals;
  218. }
  219. std::vector<std::string> readinpfile2(const std::string& fname) {
  220. std::string buf;
  221. char dump[255];
  222. std::vector<std::string> vals;
  223. std::fstream input(fname.c_str(), std::ios::in);
  224. if (input.fail()) {
  225. std::cerr << "Input file " << fname << " failed to open\n";
  226. exit(EXIT_FAILURE);
  227. }
  228. while (input >> buf) {
  229. if (buf.substr(0,2) == "//") {
  230. input.getline(dump, 255);
  231. } else {
  232. vals.push_back( std::string(buf.c_str() ));
  233. }
  234. }
  235. return vals;
  236. }