Browse Source

log spacing, but seems to reverse kernel in plots? Strange.

master
T-bone 8 years ago
parent
commit
974598a73d
3 changed files with 26 additions and 10 deletions
  1. 12
    4
      examples/KernelV0.cpp
  2. 0
    3
      include/KernelV0.h
  3. 14
    3
      src/KernelV0.cpp

+ 12
- 4
examples/KernelV0.cpp View File

41
         Kern->PushCoil( "Coil 1", Tx1 );
41
         Kern->PushCoil( "Coil 1", Tx1 );
42
         Kern->PushCoil( "Coil 2", Tx2 );
42
         Kern->PushCoil( "Coil 2", Tx2 );
43
         Kern->SetLayeredEarthEM( earth );
43
         Kern->SetLayeredEarthEM( earth );
44
-        // std::cout << *Kern << std::endl;
45
 
44
 
46
         Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
45
         Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
47
         Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
46
         Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
48
-        Kern->SetTolerance( 1e-10 );
49
-        //Kern->SetTolerance( .55 ) ; // 1%
47
+        Kern->SetTolerance( 1e-12 );
50
 
48
 
51
         Kern->SetPulseDuration(0.020);
49
         Kern->SetPulseDuration(0.020);
52
         VectorXr I(36);
50
         VectorXr I(36);
59
              10.697057141915716, 9.64778948429609, 8.709338689612677, 7.871268012862094;
57
              10.697057141915716, 9.64778948429609, 8.709338689612677, 7.871268012862094;
60
         //Kern->SetPulseCurrent( VectorXr::LinSpaced( 1, 10, 200 )  ); // nbins, low, high
58
         //Kern->SetPulseCurrent( VectorXr::LinSpaced( 1, 10, 200 )  ); // nbins, low, high
61
         Kern->SetPulseCurrent( I ); // nbins, low, high
59
         Kern->SetPulseCurrent( I ); // nbins, low, high
62
-        Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 30, 3, 45.5 ) );
60
+
61
+        //Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 30, 3, 45.5 ) ); // nlay, low, high
62
+        //10**np.linspace(np.log10(10),np.log10(19),10)
63
+        VectorXr interfaces = VectorXr::LinSpaced(31, std::log10(2), std::log10(50)); // 30 log spaced
64
+        for (int i=0; i<interfaces.size(); ++i) {
65
+            interfaces(i) = std::pow(10, interfaces(i));
66
+        }
67
+        Kern->SetDepthLayerInterfaces( interfaces ); // nlay, low, high
63
 
68
 
64
     // We could, I suppose, take the earth model in here? For non-linear that
69
     // We could, I suppose, take the earth model in here? For non-linear that
65
     // may be more natural to work with?
70
     // may be more natural to work with?
67
     std::vector<std::string> rx = {std::string("Coil 1")};
72
     std::vector<std::string> rx = {std::string("Coil 1")};
68
     Kern->CalculateK0( tx, rx, true );
73
     Kern->CalculateK0( tx, rx, true );
69
 
74
 
75
+    ofstream out = ofstream("k.yaml");
76
+    out << *Kern;
77
+    out.close();
70
 }
78
 }
71
 
79
 
72
 std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {
80
 std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {

+ 0
- 3
include/KernelV0.h View File

270
         Real                                      tol=1e-11;
270
         Real                                      tol=1e-11;
271
         Real                                      Temperature=283.;
271
         Real                                      Temperature=283.;
272
         Real                                      Taup = .020;  // Sec
272
         Real                                      Taup = .020;  // Sec
273
-        Real                                      Ip;           // Amps pulse current, deprecated PulseI
274
         Real                                      Larmor;
273
         Real                                      Larmor;
275
 
274
 
276
-        Complex                                   SUM;          // Depreceated, use Kern instead
277
-
278
         Vector3r                                  Size;
275
         Vector3r                                  Size;
279
         Vector3r                                  Origin;
276
         Vector3r                                  Origin;
280
 
277
 

+ 14
- 3
src/KernelV0.cpp View File

80
         for ( auto txm : TxRx) {
80
         for ( auto txm : TxRx) {
81
             node[txm.first] = txm.second->Serialize();
81
             node[txm.first] = txm.second->Serialize();
82
         }
82
         }
83
-
84
         // LayeredEarthEM
83
         // LayeredEarthEM
85
         node["SigmaModel"] = SigmaModel->Serialize();
84
         node["SigmaModel"] = SigmaModel->Serialize();
86
 
85
 
86
+        node["Larmor"] = Larmor;
87
+        node["Temperature"] = Temperature;
88
+        node["tol"] = tol;
89
+        node["minLevel"] = minLevel;
90
+        node["maxLevel"] = maxLevel;
91
+        node["Taup"] = Taup;
92
+
93
+        node["PulseI"] = PulseI;
94
+        node["Interfaces"] = Interfaces;
95
+
96
+        for ( int ilay=0; ilay<Interfaces.size()-1; ++ilay ) {
97
+            node["Kern-" + to_string(ilay) ] = static_cast<VectorXcr>(Kern.row(ilay));
98
+        }
87
         return node;
99
         return node;
88
     }		// -----  end of method KernelV0::Serialize  -----
100
     }		// -----  end of method KernelV0::Serialize  -----
89
 
101
 
142
         std::cout << "Calculating K0 kernel\n";
154
         std::cout << "Calculating K0 kernel\n";
143
         Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
155
         Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
144
         for (ilay=0; ilay<Interfaces.size()-1; ++ilay) {
156
         for (ilay=0; ilay<Interfaces.size()-1; ++ilay) {
145
-            std::cout << "Layer " << ilay << std::endl; //<< " q " << iq << std::endl;
157
+            std::cout << "Layer " << ilay << "\tfrom " << Interfaces(ilay) <<" to "<< Interfaces(ilay+1) << std::endl; //<< " q " << iq << std::endl;
146
             Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
158
             Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
147
             Origin(2) = Interfaces(ilay);
159
             Origin(2) = Interfaces(ilay);
148
             IntegrateOnOctreeGrid( vtkOutput );
160
             IntegrateOnOctreeGrid( vtkOutput );
162
 
174
 
163
         Vector3r cpos = Origin + Size/2.;
175
         Vector3r cpos = Origin + Size/2.;
164
 
176
 
165
-        SUM = 0;
166
         VOLSUM = 0;
177
         VOLSUM = 0;
167
         nleaves = 0;
178
         nleaves = 0;
168
         if (!vtkOutput) {
179
         if (!vtkOutput) {

Loading…
Cancel
Save