Quellcode durchsuchen

Eantenna.

pull/2/head
Trevor Irons vor 5 Jahren
Ursprung
Commit
133a27bf54
1 geänderte Dateien mit 262 neuen und 0 gelöschten Zeilen
  1. 262
    0
      Modules/FDEM1D/examples/Eantenna.cpp

+ 262
- 0
Modules/FDEM1D/examples/Eantenna.cpp Datei anzeigen

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

Laden…
Abbrechen
Speichern