mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-23 13:10:51 +08:00
Updated notebooks and minor bug fixes
This commit is contained in:
GudniRos
@@ -121,7 +121,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -132,7 +132,7 @@
|
||||
"# Define a counter\n",
|
||||
"C = simpeg.Utils.Counter()\n",
|
||||
"# Set the optimization\n",
|
||||
"opt = simpeg.Optimization.InexactGaussNewton(maxIter = 50)\n",
|
||||
"opt = simpeg.Optimization.InexactGaussNewton(maxIter = 2)\n",
|
||||
"opt.counter = C\n",
|
||||
"opt.LSshorten = 0.5\n",
|
||||
"opt.remember('xc')\n",
|
||||
@@ -155,19 +155,18 @@
|
||||
"beta = simpeg.Directives.BetaSchedule()\n",
|
||||
"betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
"targmis = simpeg.Directives.TargetMisfit()\n",
|
||||
"targmis.target = 1/2 * survey.nD\n",
|
||||
"targmis.target = survey.nD\n",
|
||||
"saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
"saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'\n",
|
||||
"# Create an inversion object\n",
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis,saveModel]) \n"
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis]) \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": false
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -178,268 +177,28 @@
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_TargMisEqnD_smoothTrue.npy'\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 3.72e+05 2.18e+05 0.00e+00 2.18e+05 4.01e+04 0 \n",
|
||||
" 1 3.72e+05 2.50e+04 2.19e-03 2.58e+04 5.62e+03 0 \n",
|
||||
" 2 3.72e+05 3.39e+03 4.60e-03 5.10e+03 1.01e+03 0 Skip BFGS \n",
|
||||
" 3 4.65e+04 1.91e+03 4.23e-03 2.10e+03 3.07e+02 0 Skip BFGS \n",
|
||||
" 4 4.65e+04 1.18e+03 1.15e-02 1.72e+03 1.62e+02 0 Skip BFGS \n",
|
||||
" 5 4.65e+04 1.00e+03 1.30e-02 1.60e+03 8.35e+01 0 \n",
|
||||
" 6 5.81e+03 8.78e+02 1.50e-02 9.66e+02 1.87e+02 0 Skip BFGS \n",
|
||||
" 7 5.81e+03 3.96e+02 5.56e-02 7.19e+02 2.51e+02 0 \n",
|
||||
" 8 5.81e+03 3.44e+02 3.87e-02 5.69e+02 1.35e+02 1 \n",
|
||||
" 9 7.26e+02 2.93e+02 4.41e-02 3.25e+02 9.13e+01 0 Skip BFGS \n",
|
||||
" 10 7.26e+02 2.41e+02 6.62e-02 2.89e+02 1.21e+02 1 \n",
|
||||
" 11 7.26e+02 1.47e+02 1.33e-01 2.44e+02 1.94e+02 0 \n",
|
||||
" 12 9.08e+01 1.35e+02 1.15e-01 1.45e+02 7.35e+01 0 \n",
|
||||
" 13 9.08e+01 8.40e+01 2.31e-01 1.05e+02 1.27e+02 0 \n",
|
||||
" 14 9.08e+01 7.09e+01 2.78e-01 9.61e+01 7.14e+01 0 \n",
|
||||
" 15 1.13e+01 6.65e+01 2.55e-01 6.94e+01 2.89e+01 0 \n",
|
||||
" 16 1.13e+01 6.41e+01 3.38e-01 6.80e+01 5.48e+01 0 \n",
|
||||
" 17 1.13e+01 5.44e+01 3.65e-01 5.85e+01 2.19e+01 0 \n",
|
||||
" 18 1.42e+00 5.04e+01 4.33e-01 5.10e+01 7.08e+01 0 Skip BFGS \n",
|
||||
" 19 1.42e+00 4.76e+01 4.52e-01 4.82e+01 1.68e+01 0 \n",
|
||||
" 20 1.42e+00 4.62e+01 4.72e-01 4.69e+01 1.39e+01 0 \n",
|
||||
" 21 1.77e-01 4.61e+01 4.68e-01 4.62e+01 1.08e+01 0 \n",
|
||||
" 22 1.77e-01 4.54e+01 4.44e-01 4.55e+01 3.34e+01 0 Skip BFGS \n",
|
||||
" 23 1.77e-01 4.40e+01 4.46e-01 4.41e+01 1.57e+01 0 \n",
|
||||
" 24 2.22e-02 4.20e+01 3.77e-01 4.20e+01 2.55e+01 2 \n",
|
||||
" 25 2.22e-02 4.03e+01 3.16e-01 4.03e+01 2.54e+01 0 \n",
|
||||
" 26 2.22e-02 3.99e+01 3.82e-01 3.99e+01 2.41e+01 0 \n",
|
||||
" 27 2.77e-03 3.93e+01 3.49e-01 3.93e+01 3.42e+01 3 Skip BFGS \n",
|
||||
" 28 2.77e-03 3.88e+01 4.32e-01 3.88e+01 2.44e+01 0 \n",
|
||||
" 29 2.77e-03 3.87e+01 4.49e-01 3.87e+01 3.02e+01 2 Skip BFGS \n",
|
||||
" 30 3.46e-04 3.86e+01 4.63e-01 3.86e+01 2.54e+01 0 \n",
|
||||
" 31 3.46e-04 3.86e+01 4.39e-01 3.86e+01 2.32e+01 0 \n",
|
||||
" 32 3.46e-04 3.84e+01 4.42e-01 3.84e+01 2.16e+01 1 Skip BFGS \n",
|
||||
" 33 4.33e-05 3.84e+01 4.54e-01 3.84e+01 2.25e+01 0 Skip BFGS \n",
|
||||
" 34 4.33e-05 3.83e+01 4.44e-01 3.83e+01 2.15e+01 0 \n",
|
||||
" 35 4.33e-05 3.72e+01 4.57e-01 3.72e+01 1.65e+01 0 \n",
|
||||
" 36 5.41e-06 3.64e+01 4.85e-01 3.64e+01 2.32e+01 1 Skip BFGS \n",
|
||||
" 37 5.41e-06 3.61e+01 4.70e-01 3.61e+01 2.34e+01 2 \n",
|
||||
" 38 5.41e-06 3.56e+01 4.35e-01 3.56e+01 2.75e+01 2 Skip BFGS \n",
|
||||
" 39 6.76e-07 3.56e+01 4.43e-01 3.56e+01 2.21e+01 0 \n",
|
||||
" 40 6.76e-07 3.55e+01 4.44e-01 3.55e+01 2.07e+01 2 \n",
|
||||
" 41 6.76e-07 3.50e+01 4.62e-01 3.50e+01 2.21e+01 2 Skip BFGS \n",
|
||||
" 42 8.45e-08 3.49e+01 4.41e-01 3.49e+01 2.31e+01 2 \n",
|
||||
" 43 8.45e-08 3.28e+01 4.83e-01 3.28e+01 3.06e+01 0 Skip BFGS \n",
|
||||
" 44 8.45e-08 3.22e+01 5.30e-01 3.22e+01 2.44e+01 1 Skip BFGS \n",
|
||||
" 45 1.06e-08 3.19e+01 5.73e-01 3.19e+01 1.06e+01 0 Skip BFGS \n",
|
||||
" 46 1.06e-08 3.18e+01 5.27e-01 3.18e+01 1.20e+01 1 \n",
|
||||
" 47 1.06e-08 3.17e+01 4.59e-01 3.17e+01 9.10e+00 0 Skip BFGS \n",
|
||||
" 48 1.32e-09 3.17e+01 5.04e-01 3.17e+01 1.08e+01 2 \n",
|
||||
" 49 1.32e-09 3.16e+01 5.22e-01 3.16e+01 9.32e+00 0 \n",
|
||||
" 50 1.32e-09 3.13e+01 4.58e-01 3.13e+01 1.42e+01 1 \n",
|
||||
" 0 2.13e+05 2.18e+05 0.00e+00 2.18e+05 4.01e+04 0 \n",
|
||||
" 1 2.13e+05 2.50e+04 2.21e-03 2.54e+04 5.65e+03 0 \n",
|
||||
" 2 2.13e+05 3.36e+03 4.80e-03 4.38e+03 9.88e+02 0 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 3.0199e-01 <= tolF*(1+|f0|) = 2.1833e+04\n",
|
||||
"1 : |xc-x_last| = 4.2553e+00 <= tolX*(1+|x0|) = 5.1104e+00\n",
|
||||
"0 : |proj(x-g)-x| = 1.4222e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 1.4222e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 50 <= iter = 50\n",
|
||||
"1 : |fc-fOld| = 2.1054e+04 <= tolF*(1+|f0|) = 2.1833e+04\n",
|
||||
"0 : |xc-x_last| = 9.1698e+00 <= tolX*(1+|x0|) = 5.1104e+00\n",
|
||||
"0 : |proj(x-g)-x| = 9.8763e+02 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 9.8763e+02 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 2 <= iter = 2\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_TargMisEqnD_smoothTrue.npy'\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 1.70e+05 2.18e+05 0.00e+00 2.18e+05 4.01e+04 0 \n",
|
||||
" 1 1.70e+05 2.50e+04 2.22e-03 2.53e+04 5.66e+03 0 \n",
|
||||
" 2 1.70e+05 3.35e+03 4.85e-03 4.18e+03 9.85e+02 0 Skip BFGS \n",
|
||||
" 3 2.13e+04 1.64e+03 5.93e-03 1.77e+03 2.67e+02 0 Skip BFGS \n",
|
||||
" 4 2.13e+04 8.27e+02 2.45e-02 1.35e+03 2.67e+02 0 \n",
|
||||
" 5 2.13e+04 7.95e+02 1.83e-02 1.18e+03 1.28e+02 0 \n",
|
||||
" 6 2.66e+03 7.09e+02 2.17e-02 7.67e+02 1.49e+02 0 \n",
|
||||
" 7 2.66e+03 4.57e+02 7.33e-02 6.52e+02 2.10e+02 0 \n",
|
||||
" 8 2.66e+03 4.28e+02 5.88e-02 5.84e+02 1.75e+02 1 \n",
|
||||
" 9 3.32e+02 3.56e+02 7.66e-02 3.82e+02 1.12e+02 0 \n",
|
||||
" 10 3.32e+02 2.71e+02 1.23e-01 3.12e+02 1.62e+02 1 \n",
|
||||
" 11 3.32e+02 2.06e+02 1.90e-01 2.69e+02 1.31e+02 1 \n",
|
||||
" 12 4.15e+01 1.69e+02 2.30e-01 1.79e+02 1.30e+02 1 \n",
|
||||
" 13 4.15e+01 1.53e+02 6.02e-01 1.78e+02 2.84e+02 0 Skip BFGS \n",
|
||||
" 14 4.15e+01 1.08e+02 4.27e-01 1.26e+02 1.01e+02 0 \n",
|
||||
" 15 5.19e+00 7.04e+01 5.63e-01 7.33e+01 7.23e+01 0 \n",
|
||||
" 16 5.19e+00 6.28e+01 6.89e-01 6.64e+01 1.31e+02 1 Skip BFGS \n",
|
||||
" 17 5.19e+00 5.21e+01 9.36e-01 5.70e+01 1.20e+02 0 Skip BFGS \n",
|
||||
" 18 6.49e-01 4.37e+01 9.80e-01 4.43e+01 7.78e+00 0 Skip BFGS \n",
|
||||
" 19 6.49e-01 4.28e+01 1.10e+00 4.35e+01 1.60e+01 0 \n",
|
||||
" 20 6.49e-01 4.27e+01 1.18e+00 4.35e+01 1.31e+01 0 \n",
|
||||
" 21 8.11e-02 4.27e+01 1.23e+00 4.28e+01 1.39e+01 2 Skip BFGS \n",
|
||||
" 22 8.11e-02 4.26e+01 1.20e+00 4.27e+01 1.26e+01 2 \n",
|
||||
" 23 8.11e-02 4.16e+01 1.47e+00 4.17e+01 1.59e+01 2 \n",
|
||||
" 24 1.01e-02 4.13e+01 1.38e+00 4.13e+01 3.20e+01 0 Skip BFGS \n",
|
||||
" 25 1.01e-02 4.10e+01 1.37e+00 4.10e+01 1.70e+01 0 \n",
|
||||
" 26 1.01e-02 4.04e+01 1.44e+00 4.05e+01 2.21e+01 1 Skip BFGS \n",
|
||||
" 27 1.27e-03 4.01e+01 1.33e+00 4.01e+01 2.03e+01 2 \n",
|
||||
" 28 1.27e-03 4.00e+01 1.49e+00 4.00e+01 2.19e+01 0 \n",
|
||||
" 29 1.27e-03 4.00e+01 1.48e+00 4.00e+01 2.50e+01 2 \n",
|
||||
" 30 1.58e-04 4.00e+01 1.49e+00 4.00e+01 2.38e+01 2 \n",
|
||||
" 31 1.58e-04 3.99e+01 1.44e+00 3.99e+01 2.75e+01 2 \n",
|
||||
" 32 1.58e-04 3.97e+01 1.64e+00 3.97e+01 2.91e+01 1 \n",
|
||||
" 33 1.98e-05 3.93e+01 1.79e+00 3.93e+01 3.19e+01 3 Skip BFGS \n",
|
||||
" 34 1.98e-05 3.90e+01 2.03e+00 3.90e+01 3.42e+01 2 Skip BFGS \n",
|
||||
" 35 1.98e-05 3.88e+01 2.44e+00 3.88e+01 3.22e+01 1 Skip BFGS \n",
|
||||
" 36 2.47e-06 3.86e+01 1.90e+00 3.86e+01 1.92e+01 0 \n",
|
||||
" 37 2.47e-06 3.85e+01 2.26e+00 3.85e+01 2.25e+01 2 \n",
|
||||
" 38 2.47e-06 3.84e+01 2.25e+00 3.84e+01 2.03e+01 0 Skip BFGS \n",
|
||||
" 39 3.09e-07 3.84e+01 2.36e+00 3.84e+01 2.10e+01 0 \n",
|
||||
" 40 3.09e-07 3.84e+01 2.49e+00 3.84e+01 2.08e+01 2 Skip BFGS \n",
|
||||
" 41 3.09e-07 3.84e+01 2.37e+00 3.84e+01 2.08e+01 1 \n",
|
||||
" 42 3.87e-08 3.78e+01 2.59e+00 3.78e+01 1.62e+01 1 Skip BFGS \n",
|
||||
" 43 3.87e-08 3.77e+01 2.62e+00 3.77e+01 1.14e+01 0 Skip BFGS \n",
|
||||
" 44 3.87e-08 3.76e+01 2.83e+00 3.76e+01 1.67e+01 3 \n",
|
||||
" 45 4.83e-09 3.75e+01 2.80e+00 3.75e+01 2.17e+01 2 Skip BFGS \n",
|
||||
" 46 4.83e-09 3.73e+01 2.73e+00 3.73e+01 1.19e+01 0 \n",
|
||||
" 47 4.83e-09 3.71e+01 2.46e+00 3.71e+01 1.67e+01 2 \n",
|
||||
" 48 6.04e-10 3.71e+01 2.33e+00 3.71e+01 1.45e+01 0 \n",
|
||||
" 49 6.04e-10 3.71e+01 2.46e+00 3.71e+01 1.85e+01 1 \n",
|
||||
" 50 6.04e-10 3.70e+01 2.39e+00 3.70e+01 1.99e+01 3 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 2.5694e-02 <= tolF*(1+|f0|) = 2.1833e+04\n",
|
||||
"1 : |xc-x_last| = 1.9839e+00 <= tolX*(1+|x0|) = 5.1104e+00\n",
|
||||
"0 : |proj(x-g)-x| = 1.9850e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 1.9850e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 50 <= iter = 50\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_TargMisEqnD_smoothTrue.npy'\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 1.19e+05 2.18e+05 0.00e+00 2.18e+05 4.01e+04 0 Skip BFGS \n",
|
||||
" 1 1.19e+05 2.50e+04 2.22e-03 2.52e+04 5.67e+03 0 Skip BFGS \n",
|
||||
" 2 1.19e+05 3.35e+03 4.92e-03 3.93e+03 9.84e+02 0 Skip BFGS \n",
|
||||
" 3 1.49e+04 1.53e+03 7.13e-03 1.63e+03 2.55e+02 0 Skip BFGS \n",
|
||||
" 4 1.49e+04 7.29e+02 3.29e-02 1.22e+03 2.94e+02 0 \n",
|
||||
" 5 1.49e+04 6.80e+02 2.16e-02 1.00e+03 1.41e+02 0 \n",
|
||||
" 6 1.86e+03 6.17e+02 2.45e-02 6.63e+02 1.44e+02 0 \n",
|
||||
" 7 1.86e+03 4.77e+02 7.62e-02 6.19e+02 2.19e+02 0 \n",
|
||||
" 8 1.86e+03 4.43e+02 6.42e-02 5.63e+02 2.03e+02 1 \n",
|
||||
" 9 2.32e+02 3.74e+02 8.67e-02 3.94e+02 1.32e+02 0 \n",
|
||||
" 10 2.32e+02 2.67e+02 1.91e-01 3.11e+02 1.89e+02 1 \n",
|
||||
" 11 2.32e+02 1.29e+02 5.16e-01 2.49e+02 1.79e+02 0 Skip BFGS \n",
|
||||
" 12 2.91e+01 1.32e+02 3.87e-01 1.43e+02 1.57e+02 1 \n",
|
||||
" 13 2.91e+01 1.11e+02 6.45e-01 1.30e+02 1.91e+02 1 \n",
|
||||
" 14 2.91e+01 6.74e+01 9.75e-01 9.57e+01 1.25e+02 0 Skip BFGS \n",
|
||||
" 15 3.63e+00 5.28e+01 1.03e+00 5.66e+01 1.14e+01 0 Skip BFGS \n",
|
||||
" 16 3.63e+00 4.93e+01 1.04e+00 5.31e+01 3.75e+01 1 Skip BFGS \n",
|
||||
" 17 3.63e+00 4.65e+01 1.16e+00 5.07e+01 8.86e+01 0 \n",
|
||||
" 18 4.54e-01 4.36e+01 1.16e+00 4.42e+01 5.66e+01 0 \n",
|
||||
" 19 4.54e-01 4.18e+01 1.21e+00 4.23e+01 5.22e+01 1 Skip BFGS \n",
|
||||
" 20 4.54e-01 3.95e+01 1.24e+00 4.01e+01 1.23e+01 0 \n",
|
||||
" 21 5.68e-02 3.93e+01 1.26e+00 3.94e+01 1.03e+01 0 \n",
|
||||
" 22 5.68e-02 3.92e+01 1.29e+00 3.93e+01 9.85e+00 0 \n",
|
||||
" 23 5.68e-02 3.89e+01 1.44e+00 3.90e+01 1.92e+01 1 Skip BFGS \n",
|
||||
" 24 7.09e-03 3.87e+01 1.36e+00 3.87e+01 1.35e+01 0 \n",
|
||||
" 25 7.09e-03 3.86e+01 1.34e+00 3.86e+01 1.28e+01 0 Skip BFGS \n",
|
||||
" 26 7.09e-03 3.85e+01 1.36e+00 3.85e+01 1.29e+01 0 \n",
|
||||
" 27 8.87e-04 3.79e+01 1.32e+00 3.79e+01 1.64e+01 1 Skip BFGS \n",
|
||||
" 28 8.87e-04 3.79e+01 1.30e+00 3.79e+01 1.60e+01 0 \n",
|
||||
" 29 8.87e-04 3.79e+01 1.22e+00 3.79e+01 2.12e+01 0 Skip BFGS \n",
|
||||
" 30 1.11e-04 3.77e+01 1.26e+00 3.77e+01 1.51e+01 0 \n",
|
||||
" 31 1.11e-04 3.77e+01 1.24e+00 3.77e+01 1.43e+01 0 Skip BFGS \n",
|
||||
" 32 1.11e-04 3.76e+01 1.24e+00 3.76e+01 1.51e+01 0 \n",
|
||||
" 33 1.39e-05 3.75e+01 1.29e+00 3.75e+01 1.85e+01 2 Skip BFGS \n",
|
||||
" 34 1.39e-05 3.75e+01 1.24e+00 3.75e+01 1.58e+01 1 \n",
|
||||
" 35 1.39e-05 3.75e+01 1.17e+00 3.75e+01 1.76e+01 1 Skip BFGS \n",
|
||||
" 36 1.73e-06 3.75e+01 1.24e+00 3.75e+01 1.61e+01 1 \n",
|
||||
" 37 1.73e-06 3.74e+01 1.41e+00 3.74e+01 1.76e+01 3 Skip BFGS \n",
|
||||
" 38 1.73e-06 3.74e+01 1.37e+00 3.74e+01 1.76e+01 2 \n",
|
||||
" 39 2.17e-07 3.73e+01 1.45e+00 3.73e+01 2.92e+01 1 Skip BFGS \n",
|
||||
" 40 2.17e-07 3.71e+01 1.51e+00 3.71e+01 2.50e+01 1 \n",
|
||||
" 41 2.17e-07 3.70e+01 1.46e+00 3.70e+01 3.25e+01 0 \n",
|
||||
" 42 2.71e-08 3.67e+01 1.40e+00 3.67e+01 3.70e+01 2 Skip BFGS \n",
|
||||
" 43 2.71e-08 3.67e+01 1.39e+00 3.67e+01 3.64e+01 2 Skip BFGS \n",
|
||||
" 44 2.71e-08 3.66e+01 1.40e+00 3.66e+01 3.47e+01 2 \n",
|
||||
" 45 3.38e-09 3.62e+01 1.30e+00 3.62e+01 2.46e+01 1 Skip BFGS \n",
|
||||
" 46 3.38e-09 3.61e+01 1.24e+00 3.61e+01 1.97e+01 1 Skip BFGS \n",
|
||||
" 47 3.38e-09 3.61e+01 1.28e+00 3.61e+01 2.11e+01 1 \n",
|
||||
" 48 4.23e-10 3.60e+01 1.28e+00 3.60e+01 1.81e+01 0 Skip BFGS \n",
|
||||
" 49 4.23e-10 3.60e+01 1.26e+00 3.60e+01 1.81e+01 1 \n",
|
||||
" 50 4.23e-10 3.59e+01 1.27e+00 3.59e+01 1.80e+01 1 \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 1.8765e-02 <= tolF*(1+|f0|) = 2.1833e+04\n",
|
||||
"1 : |xc-x_last| = 4.3809e-01 <= tolX*(1+|x0|) = 5.1104e+00\n",
|
||||
"0 : |proj(x-g)-x| = 1.8033e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 1.8033e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 50 <= iter = 50\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_TargMisEqnD_smoothTrue.npy'\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 2.04e+05 2.18e+05 0.00e+00 2.18e+05 4.01e+04 0 \n",
|
||||
" 1 2.04e+05 2.50e+04 2.21e-03 2.54e+04 5.65e+03 0 \n",
|
||||
" 2 2.04e+05 3.36e+03 4.81e-03 4.34e+03 9.87e+02 0 Skip BFGS \n",
|
||||
" 3 2.55e+04 1.70e+03 5.43e-03 1.84e+03 2.74e+02 0 Skip BFGS \n",
|
||||
" 4 2.55e+04 8.94e+02 2.07e-02 1.42e+03 2.45e+02 0 \n",
|
||||
" 5 2.55e+04 8.27e+02 1.76e-02 1.28e+03 1.27e+02 0 \n",
|
||||
" 6 3.19e+03 7.66e+02 1.96e-02 8.29e+02 1.52e+02 0 \n",
|
||||
" 7 3.19e+03 4.82e+02 7.40e-02 7.18e+02 2.52e+02 0 \n",
|
||||
" 8 3.19e+03 4.13e+02 6.65e-02 6.25e+02 1.58e+02 0 \n",
|
||||
" 9 3.98e+02 3.81e+02 7.32e-02 4.11e+02 1.43e+02 0 \n",
|
||||
" 10 3.98e+02 2.69e+02 1.28e-01 3.20e+02 1.49e+02 1 \n",
|
||||
" 11 3.98e+02 2.01e+02 2.10e-01 2.84e+02 1.84e+02 0 \n",
|
||||
" 12 4.98e+01 1.66e+02 2.25e-01 1.78e+02 1.25e+02 0 \n",
|
||||
" 13 4.98e+01 1.28e+02 3.65e-01 1.46e+02 9.89e+01 0 \n",
|
||||
" 14 4.98e+01 1.05e+02 4.50e-01 1.27e+02 1.64e+02 1 \n",
|
||||
" 15 6.23e+00 6.86e+01 5.39e-01 7.19e+01 1.64e+02 0 Skip BFGS \n",
|
||||
" 16 6.23e+00 4.97e+01 5.96e-01 5.34e+01 1.16e+01 0 Skip BFGS \n",
|
||||
" 17 6.23e+00 4.83e+01 6.74e-01 5.25e+01 1.77e+01 0 Skip BFGS \n",
|
||||
" 18 7.78e-01 4.76e+01 6.83e-01 4.81e+01 8.50e+00 0 \n",
|
||||
" 19 7.78e-01 4.67e+01 7.68e-01 4.73e+01 6.91e+00 0 Skip BFGS \n",
|
||||
" 20 7.78e-01 4.65e+01 7.23e-01 4.71e+01 3.93e+01 0 \n",
|
||||
" 21 9.73e-02 4.58e+01 8.12e-01 4.59e+01 5.65e+01 2 \n",
|
||||
" 22 9.73e-02 4.55e+01 7.62e-01 4.55e+01 3.81e+01 0 \n",
|
||||
" 23 9.73e-02 4.52e+01 6.89e-01 4.53e+01 3.60e+01 2 \n",
|
||||
" 24 1.22e-02 4.48e+01 7.67e-01 4.48e+01 3.68e+01 1 \n",
|
||||
" 25 1.22e-02 4.46e+01 7.74e-01 4.46e+01 3.70e+01 3 \n",
|
||||
" 26 1.22e-02 4.45e+01 7.45e-01 4.45e+01 4.29e+01 2 \n",
|
||||
" 27 1.52e-03 4.43e+01 7.79e-01 4.43e+01 4.13e+01 2 \n",
|
||||
" 28 1.52e-03 4.39e+01 7.69e-01 4.39e+01 4.83e+01 2 \n",
|
||||
" 29 1.52e-03 4.38e+01 7.82e-01 4.38e+01 6.21e+01 0 Skip BFGS \n",
|
||||
" 30 1.90e-04 4.36e+01 7.74e-01 4.36e+01 5.22e+01 1 \n",
|
||||
" 31 1.90e-04 4.36e+01 7.93e-01 4.36e+01 5.35e+01 2 \n",
|
||||
" 32 1.90e-04 4.31e+01 7.77e-01 4.31e+01 5.43e+01 1 \n",
|
||||
" 33 2.38e-05 4.11e+01 8.44e-01 4.11e+01 2.16e+01 0 \n",
|
||||
" 34 2.38e-05 4.05e+01 9.56e-01 4.05e+01 2.25e+01 0 \n",
|
||||
" 35 2.38e-05 4.03e+01 9.07e-01 4.03e+01 1.80e+01 1 \n",
|
||||
" 36 2.97e-06 4.03e+01 8.80e-01 4.03e+01 2.09e+01 3 \n",
|
||||
" 37 2.97e-06 4.03e+01 8.13e-01 4.03e+01 2.11e+01 2 Skip BFGS \n",
|
||||
" 38 2.97e-06 4.01e+01 8.53e-01 4.01e+01 2.04e+01 2 \n",
|
||||
" 39 3.71e-07 4.00e+01 8.76e-01 4.00e+01 1.56e+01 0 \n",
|
||||
" 40 3.71e-07 3.99e+01 9.15e-01 3.99e+01 1.25e+01 0 \n",
|
||||
" 41 3.71e-07 3.94e+01 9.05e-01 3.94e+01 1.02e+01 1 \n",
|
||||
" 42 4.64e-08 3.92e+01 9.48e-01 3.92e+01 8.52e+00 0 Skip BFGS \n",
|
||||
" 43 4.64e-08 3.91e+01 8.92e-01 3.91e+01 4.37e+00 0 \n",
|
||||
" 44 4.64e-08 3.89e+01 9.26e-01 3.89e+01 1.63e+01 1 \n",
|
||||
" 45 5.80e-09 3.88e+01 8.77e-01 3.88e+01 1.21e+01 0 \n",
|
||||
" 46 5.80e-09 3.88e+01 9.08e-01 3.88e+01 1.54e+01 0 \n",
|
||||
" 47 5.80e-09 3.83e+01 7.94e-01 3.83e+01 1.77e+01 0 \n",
|
||||
" 48 7.25e-10 3.81e+01 8.11e-01 3.81e+01 1.68e+01 1 \n",
|
||||
" 49 7.25e-10 3.81e+01 7.94e-01 3.81e+01 1.47e+01 0 \n",
|
||||
" 50 7.25e-10 3.80e+01 8.16e-01 3.80e+01 1.90e+01 3 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 1.5226e-01 <= tolF*(1+|f0|) = 2.1833e+04\n",
|
||||
"1 : |xc-x_last| = 1.1100e+00 <= tolX*(1+|x0|) = 5.1104e+00\n",
|
||||
"0 : |proj(x-g)-x| = 1.8975e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 1.8975e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 50 <= iter = 50\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"1 loops, best of 3: 20min 52s per loop\n"
|
||||
" "
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%%timeit\n",
|
||||
"##### \n",
|
||||
"# Run the inversion, given the background model as a start.\n",
|
||||
"mopt = inv.run(m_0)"
|
||||
"import cProfile\n",
|
||||
"%prun mopt = inv.run(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
File diff suppressed because one or more lines are too long
@@ -21,6 +21,7 @@
|
||||
"source": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n",
|
||||
"from simpegMT.Utils.dataUtils import rec2ndarr\n",
|
||||
"import cPickle as pickle"
|
||||
]
|
||||
},
|
||||
@@ -62,6 +63,61 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"995.625"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"(np.sum(bPad)+600)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ 151.875, 101.25 , 67.5 , 45. , 30. ])"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"bPad"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"!nautilus ."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Load the model to the uniform cell mesh\n",
|
||||
@@ -72,7 +128,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
@@ -87,7 +143,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -101,26 +157,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "SolverException",
|
||||
"evalue": "Mumps Exception [-13] - An error occurred in a Fortran ALLOCATE statement. The size that the package requested is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million.",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[1;31mSolverException\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[1;32m<ipython-input-7-10aaec78a1f3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 21\u001b[0m \u001b[1;31m# Forward model the data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 22\u001b[1;33m \u001b[0mfields\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodelTD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 23\u001b[0m \u001b[0mmtData\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfields\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT3D/Problems.pyc\u001b[0m in \u001b[0;36mfields\u001b[1;34m(self, m)\u001b[0m\n\u001b[0;32m 100\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetA\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 101\u001b[0m \u001b[0mrhs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetRHS\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 102\u001b[1;33m \u001b[0mAinv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSolver\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msolverOpts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 103\u001b[0m \u001b[0me_s\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mAinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mrhs\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 104\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/pymatsolver/pymatsolver/Mumps/__init__.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, A, symmetric, fromPointer)\u001b[0m\n\u001b[0;32m 96\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 97\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mfromPointer\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 98\u001b[1;33m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfactor\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 99\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfromPointer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_Pointer\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 100\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpointer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfromPointer\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/pymatsolver/pymatsolver/Mumps/__init__.pyc\u001b[0m in \u001b[0;36mfactor\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 132\u001b[0m self.A.indptr+1)\n\u001b[0;32m 133\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mierr\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 134\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mSolverException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Mumps Exception [%d] - %s\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_mumpsErrors\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 135\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0mierr\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 136\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"Mumps Warning [%d] - %s\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_mumpsErrors\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;31mSolverException\u001b[0m: Mumps Exception [-13] - An error occurred in a Fortran ALLOCATE statement. The size that the package requested is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million."
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Make the receiver list\n",
|
||||
"rxList = []\n",
|
||||
@@ -134,50 +175,340 @@
|
||||
"# Survey MT\n",
|
||||
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
|
||||
"\n",
|
||||
"# Setup the problem object\n",
|
||||
"sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model\n",
|
||||
"problem = simpegmt.ProblemMT3D.eForm_ps(mesh3d,sigmaPrimary = sigma1d)\n",
|
||||
"problem.verbose = False\n",
|
||||
"from pymatsolver import MumpsSolver\n",
|
||||
"problem.Solver = MumpsSolver\n",
|
||||
"problem.pair(survey)\n",
|
||||
"\n",
|
||||
"# Forward model the data\n",
|
||||
"fields = problem.fields(modelTD)\n",
|
||||
"mtData = survey.projectFields(fields)"
|
||||
"if False:\n",
|
||||
" fields = problem.fields(modelTD)\n",
|
||||
" mtData = survey.projectFields(fields)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2834: FutureWarning: Numpy has detected that you (may be) writing to an array returned\n",
|
||||
"by numpy.diagonal or by selecting multiple fields in a record\n",
|
||||
"array. This code will likely break in a future numpy release --\n",
|
||||
"see numpy.diagonal or arrays.indexing reference docs for details.\n",
|
||||
"The quick fix is to make an explicit copy (e.g., do\n",
|
||||
"arr.diagonal().copy() or arr[['f0','f1']].copy()).\n",
|
||||
" if (obj.__array_interface__[\"data\"][0]\n",
|
||||
"/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2835: FutureWarning: Numpy has detected that you (may be) writing to an array returned\n",
|
||||
"by numpy.diagonal or by selecting multiple fields in a record\n",
|
||||
"array. This code will likely break in a future numpy release --\n",
|
||||
"see numpy.diagonal or arrays.indexing reference docs for details.\n",
|
||||
"The quick fix is to make an explicit copy (e.g., do\n",
|
||||
"arr.diagonal().copy() or arr[['f0','f1']].copy()).\n",
|
||||
" != self.__array_interface__[\"data\"][0]):\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Load the data\n",
|
||||
"mtSeogiData = simpegmt.DataMT.DataMT(survey,np.load('seogiModel_MTdata.npy'))\n",
|
||||
"mtSeogirecData = mtSeogiData.toRecArray('Complex')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
" 50*50*48"
|
||||
"sys.path.append('/home/gudni/Dropbox/code/python/MTview/')\n",
|
||||
"import interactivePlotFunctions as iPf"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"np.sum(modTD<1e-7)"
|
||||
"% matplotlib qt\n",
|
||||
"# Plot the data\n",
|
||||
"mtSeogiMap = iPf.MTinteractiveMap([mtSeogirecData])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"45000/50**2"
|
||||
"# Add noise to the data\n",
|
||||
"std = 0.05 # 5% std\n",
|
||||
"if os.path.isfile('seogiModel_MTdata.npy') and os.path.isfile('seogiModel_dobs.npy'):\n",
|
||||
" d_true = np.load('seogiModel_MTdata.npy')\n",
|
||||
" d_obs = np.load('seogiModel_dobs.npy')\n",
|
||||
"else:\n",
|
||||
" d_true = simpeg.mkvc(mtSeogiData)\n",
|
||||
" d_obs = d_true + std*abs(d_true)*np.random.randn(*d_true.shape)\n",
|
||||
" np.save('seogiModel_dobs.npy',d_obs)\n",
|
||||
"# Assign the dobs\n",
|
||||
"survey.dtrue = d_true\n",
|
||||
"survey.dobs = d_obs\n",
|
||||
"survey.std = np.abs(survey.dobs*std) + 0.01*np.linalg.norm(survey.dobs) #survey.dobs*0 + std\n",
|
||||
"# Assign the data weight\n",
|
||||
"survey.Wd = 1/survey.std #(abs(survey.dobs)*survey.std)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Make an MTdata object with dobs\n",
|
||||
"mtSeogiDobs = simpegmt.DataMT.DataMT(survey,survey.dobs)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 36,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Run the setup\n",
|
||||
"mt1DdataZyxList = simpegmt.Utils.dataUtils.convert3Dto1Dobject(mtSeogiDobs,'zyx')\n",
|
||||
"mt1DdataZxyList = simpegmt.Utils.dataUtils.convert3Dto1Dobject(mtSeogiDobs,'zxy')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"m1d = simpeg.Mesh.TensorMesh([mesh3d.hz],[mesh3d.x0[-1]])\n",
|
||||
"# Setup the problem object\n",
|
||||
"sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model\n",
|
||||
"# Set the mapping\n",
|
||||
"actMap = simpeg.Maps.ActiveCells(m1d, sigma1d > 1e-7, np.log(1e-8), nC=m1d.nCx)\n",
|
||||
"mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap\n",
|
||||
"problem = simpegmt.ProblemMT1D.eForm_psField(m1d,sigmaPrimary = sigma1d,mapping=mappingExpAct)\n",
|
||||
"problem.verbose = False\n",
|
||||
"from pymatsolver import MumpsSolver\n",
|
||||
"problem.Solver = MumpsSolver\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([[-1989.0625 , -1989.0625 , -995.625 ],\n",
|
||||
" [-1609.375 , -1989.0625 , -995.625 ],\n",
|
||||
" [-1356.25 , -1989.0625 , -995.625 ],\n",
|
||||
" ..., \n",
|
||||
" [ 1356.25 , 1989.0625 , 2538.2509238],\n",
|
||||
" [ 1609.375 , 1989.0625 , 2538.2509238],\n",
|
||||
" [ 1989.0625 , 1989.0625 , 2538.2509238]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"mesh3d.gridN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 48,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Define a function to run the forward problem\n",
|
||||
"def runInversionModel(data,problem,m1d,nameflag):\n",
|
||||
" problem.unpair()\n",
|
||||
" data.survey.unpair()\n",
|
||||
" problem.pair(data.survey)\n",
|
||||
" # Define a counter\n",
|
||||
" C = simpeg.Utils.Counter()\n",
|
||||
" # Set the optimization\n",
|
||||
" opt = simpeg.Optimization.InexactGaussNewton(maxIter = 20)\n",
|
||||
" opt.counter = C\n",
|
||||
" opt.LSshorten = 0.5\n",
|
||||
" opt.remember('xc')\n",
|
||||
" # Data misfit\n",
|
||||
" dmis = simpeg.DataMisfit.l2_DataMisfit(data.survey)\n",
|
||||
" # Regularization\n",
|
||||
" regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)\n",
|
||||
" reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
" reg.smoothModel = False\n",
|
||||
" reg.alpha_s = 1e-6\n",
|
||||
" reg.alpha_x = 1.\n",
|
||||
" # reg.alpha_xx = .001\n",
|
||||
" # Inversion problem\n",
|
||||
" invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)\n",
|
||||
" invProb.counter = C\n",
|
||||
" # Beta cooling\n",
|
||||
" beta = simpeg.Directives.BetaSchedule()\n",
|
||||
" betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
" targmis = simpeg.Directives.TargetMisfit()\n",
|
||||
" targmis.target = .6 * data.survey.nD\n",
|
||||
" print 'Target misfit is {:.0f}'.format(targmis.target)\n",
|
||||
" locs = np.ones((regMesh.nC,1))*data.survey.srcList[0].rxList[0].locs[:,:-1]\n",
|
||||
" saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
" saveModel.fileName = 'Inversion_modAt{:.0f}_{:.0f}'.format(locs[0,0],locs[0,1])\n",
|
||||
" print saveModel.fileName\n",
|
||||
" # Create an inversion object\n",
|
||||
" inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])\n",
|
||||
" # Run\n",
|
||||
" m_0 = np.log(1e-2+0*problem.sigmaPrimary[problem.mapping.sigmaMap.maps[-1].indActive])\n",
|
||||
" problem.survey.mtrue = m_0\n",
|
||||
" mopt = inv.run(m_0)\n",
|
||||
" # Save the model\n",
|
||||
" \n",
|
||||
" xyzv = np.hstack((locs,simpeg.mkvc(regMesh.gridCC,2),simpeg.mkvc(mopt,2)))\n",
|
||||
" np.save('xyzmod_{:s}_{:.0f}_{:.0f}.npy'.format(nameflag,locs[0,0],locs[0,1]),xyzv)\n",
|
||||
"# return mopt\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 42,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# print np.exp(problem.survey.mtrue)\n",
|
||||
"# print problem.survey.nD"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 49,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Target misfit is 16\n",
|
||||
"SimPEG.InvProblem will set Regularization.mref to m0.\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 1.65e+03 9.20e+02 0.00e+00 9.20e+02 4.30e+02 0 \n",
|
||||
" 1 1.65e+03 5.03e+02 6.02e-04 5.04e+02 6.99e+02 0 \n",
|
||||
" 2 1.65e+03 3.24e+01 3.37e-04 3.29e+01 7.11e+01 0 \n",
|
||||
" 3 2.06e+02 1.72e+01 1.75e-03 1.76e+01 2.15e+01 0 Skip BFGS \n",
|
||||
" 4 2.06e+02 1.64e+01 2.81e-03 1.70e+01 2.36e+01 1 \n",
|
||||
" 5 2.06e+02 1.62e+01 3.17e-03 1.69e+01 2.30e+01 1 \n",
|
||||
" 6 2.58e+01 1.62e+01 3.13e-03 1.63e+01 2.27e+01 1 \n",
|
||||
" 7 2.58e+01 1.61e+01 5.41e-03 1.63e+01 2.77e+01 2 \n",
|
||||
" 8 2.58e+01 1.61e+01 6.81e-03 1.63e+01 3.16e+01 2 \n",
|
||||
" 9 3.23e+00 1.60e+01 6.54e-03 1.60e+01 3.18e+01 3 \n",
|
||||
" 10 3.23e+00 1.57e+01 1.13e-02 1.58e+01 3.16e+01 3 \n",
|
||||
" 11 3.23e+00 1.56e+01 8.83e-03 1.57e+01 3.22e+01 3 \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 9.2072e+01\n",
|
||||
"1 : |xc-x_last| = 1.5187e-01 <= tolX*(1+|x0|) = 2.6224e+00\n",
|
||||
"0 : |proj(x-g)-x| = 3.2196e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 3.2196e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"0 : maxIter = 20 <= iter = 12\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"for dat in mt1DdataZyxList[0:1]:\n",
|
||||
" runInversionModel(dat,problem,m1d,'zyx')\n",
|
||||
"\n",
|
||||
"# for dat in mt1DdataZxyList:\n",
|
||||
"# runInversionModel(dat,problem,m1d,'zxy')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 44,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"xyzmod_zxy_-700_-500.npy xyzmod_zyx_-700_-500.npy\r\n",
|
||||
"xyzmod_zxy_-700_-700.npy xyzmod_zyx_-700_-700.npy\r\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"ls xyzmod*"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 46,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ -919.6875 , -793.125 , -708.75 , -652.5 ,\n",
|
||||
" -615. , -590. , -570. , -550. ,\n",
|
||||
" -530. , -510. , -490. , -470. ,\n",
|
||||
" -450. , -430. , -410. , -390. ,\n",
|
||||
" -370. , -350. , -330. , -310. ,\n",
|
||||
" -290. , -270. , -250. , -230. ,\n",
|
||||
" -210. , -190. , -170. , -150. ,\n",
|
||||
" -130. , -110. , -90. , -70. ,\n",
|
||||
" -50. , -30. , -10. , 13. ,\n",
|
||||
" 42.9 , 81.77 , 132.301 , 197.9913 ,\n",
|
||||
" 283.38869 , 394.405297 , 538.7268861 , 726.34495193,\n",
|
||||
" 970.24843751, 1287.32296876, 1699.51985939, 2235.37581721])"
|
||||
]
|
||||
},
|
||||
"execution_count": 46,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"m1d.gridCC\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -187,7 +518,184 @@
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
"source": [
|
||||
"dview = rc[:]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"dview['problem'] = problem\n",
|
||||
"dview['m1d'] = m1d"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"dview.map_sync(lambda l:runInversionModel(l,problem,m1d),mt1DdataList)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"m_0 = np.log(problem.sigmaPrimary[problem.mapping.sigmaMap.maps[-1].indActive])\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,[m_0,mopt])\n",
|
||||
"plt.suptitle('Target misfit-smooth False')\n",
|
||||
"plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"np.exp(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import multiprocessing, os, numpy as np\n",
|
||||
"from glob import glob\n",
|
||||
"import subprocess\n",
|
||||
"\n",
|
||||
"def runMT3Dfwd(path):\n",
|
||||
" \"\"\" Worker function \"\"\"\n",
|
||||
" orgDir = os.getcwd()\n",
|
||||
" print 'Starting in {:s}\\n'.format(orgDir)\n",
|
||||
" \n",
|
||||
" os.chdir(path)\n",
|
||||
" print '########### \\nRunning in {:s}\\n###########'.format(os.getcwd())\n",
|
||||
" cmd = '/tera_raid/gudni/Codes/ZTEM_MT3Dinv_06082014/MT3Dfwd'\n",
|
||||
" # os.system(cmd)\n",
|
||||
" subprocess.call(cmd)\n",
|
||||
" os.chdir(orgDir)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"if __name__ == '__main__':\n",
|
||||
" pool = multiprocessing.Pool(processes=12)\n",
|
||||
" foldList = glob('MT3Dfwd/model100w*/Freq*')\n",
|
||||
" \n",
|
||||
" pool.map(runMT3Dfwd,foldList)\n",
|
||||
" pool.cself.mesh.getEdgeInnerProductlose()\n",
|
||||
" pool.join()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"locs=rec2ndarr(uniLocs)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"simpeg.mkvc(locs[0,:],2).T"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"simpegmt.DataMT.DataMT"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"sys.path.append('/home/gudni/Dropbox/code/python/MTview/')\n",
|
||||
"import interactivePlotFunctions as iPf\n",
|
||||
"dataRec = mt1DdataList[0].toRecArray('Complex')\n",
|
||||
"iPf.MTinteractiveMap([dataRec])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"simpegmt.Utils.dataUtils.getAppRes(mt1DdataList[0])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
" "
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
||||
@@ -9,9 +9,7 @@ from pymatsolver import MumpsSolver
|
||||
import simpegMT as simpegmt, SimPEG as simpeg
|
||||
import numpy as np, scipy
|
||||
## Setup the forward modeling
|
||||
# Read the model
|
||||
modelname = "simpegTDmodel.con"
|
||||
sigma = np.loadtxt(modelname)
|
||||
|
||||
# Make the mesh.
|
||||
mTensor = simpeg.Utils.meshTensor
|
||||
cSize = [50,20]
|
||||
@@ -23,20 +21,27 @@ x0 = np.array([-1250,-1250,- 30*20])
|
||||
mesh3dCons = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)
|
||||
# With padding
|
||||
hPad = mTensor([(cSize[0],5,1.5)])
|
||||
aPad = mTensor([(cSize[1],13,1.3)])
|
||||
bPad = mTensor([(cSize[1],5,-1.5)])
|
||||
aPad = mTensor([(cSize[1],9,1.5)])
|
||||
bPad = mTensor([(cSize[1],9,-1.5)])
|
||||
hxPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))
|
||||
hyPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))
|
||||
hzPad = np.hstack((bPad,mTensor([(cSize[1],30)]),aPad))
|
||||
x0Pad = np.array([-(np.sum(hPad)+1000),-(np.sum(hPad)+1000),-(np.sum(bPad)+600)])
|
||||
x0Pad = np.array([-(np.sum(hPad)+1000),-(np.sum(hPad)+1000),-(np.sum(bPad)+(20*30))])
|
||||
mesh3d = simpeg.Mesh.TensorMesh([hxPad,hyPad,hzPad],x0Pad)
|
||||
|
||||
# Read the model
|
||||
modelname = "simpegTDmodel.con"
|
||||
|
||||
# Load the model to the uniform cell mesh
|
||||
modelUniCell = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3dCons)
|
||||
|
||||
# Load the model to the mesh with padding cells
|
||||
modelTD = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)
|
||||
|
||||
modelT = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)
|
||||
# Adjust the model to reflect changes in the mesh (fewer aircells)
|
||||
modMat = mesh3d.r(modelT,'CC','CC','M')
|
||||
modNewMat = np.ones((50,50,48))*modMat[0,0,0]
|
||||
modNewMat[:,:,9::] = modMat[:,:,:-9]
|
||||
modelTD = mesh3d.r(modNewMat,'CC','CC','V')
|
||||
|
||||
# Define the data locations
|
||||
xG,yG = np.meshgrid(np.linspace(-700,700,8),np.linspace(-700,700,8))
|
||||
@@ -58,7 +63,7 @@ survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
# Setup the problem object
|
||||
sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model
|
||||
problem = simpegmt.ProblemMT3D.eForm_ps(mesh3d,sigmaPrimary = sigma1d)
|
||||
problem.verbose = False
|
||||
problem.verbose = True
|
||||
from pymatsolver import MumpsSolver
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
|
||||
+1
-1
@@ -91,7 +91,7 @@ class DataMT(Survey.Data):
|
||||
if srcType=='primary':
|
||||
src = simpegMT.SurveyMT.srcMT_polxy_1Dprimary
|
||||
elif srcType=='total':
|
||||
simpegMT.SurveyMT.srcMT_polxy_1DhomotD
|
||||
src = simpegMT.SurveyMT.srcMT_polxy_1DhomotD
|
||||
else:
|
||||
raise NotImplementedError('{:s} is not a valid source type for MTdata')
|
||||
|
||||
|
||||
+11
-10
@@ -89,8 +89,8 @@ class RxMT(Survey.BaseRx):
|
||||
'''
|
||||
|
||||
if self.projType is 'Z1D':
|
||||
Pex = mesh.getInterpolationMat(self.locs,'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
|
||||
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
|
||||
ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
# Note: Has a minus sign in front, to comply with quadrant calculations.
|
||||
@@ -144,8 +144,8 @@ class RxMT(Survey.BaseRx):
|
||||
|
||||
if not adjoint:
|
||||
if self.projType is 'Z1D':
|
||||
Pex = mesh.getInterpolationMat(self.locs,'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
|
||||
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
|
||||
@@ -160,8 +160,8 @@ class RxMT(Survey.BaseRx):
|
||||
elif adjoint:
|
||||
# Note: The v vector is real and the return should be complex
|
||||
if self.projType is 'Z1D':
|
||||
Pex = mesh.getInterpolationMat(self.locs,'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
|
||||
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
|
||||
@@ -223,14 +223,15 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
|
||||
self.sigma1d = None
|
||||
srcMT.__init__(self, rxList, freq)
|
||||
|
||||
|
||||
# Hidden property of the ePrimary
|
||||
self._ePrimary = None
|
||||
|
||||
def ePrimary(self,problem):
|
||||
# Get primary fields for both polarizations
|
||||
self.sigma1d = problem._sigmaPrimary
|
||||
eBG_bp = homo1DModelSource(problem.mesh,self.freq,self.sigma1d)
|
||||
return eBG_bp
|
||||
if self._ePrimary is None:
|
||||
self._ePrimary = homo1DModelSource(problem.mesh,self.freq,self.sigma1d)
|
||||
return self._ePrimary
|
||||
|
||||
def bPrimary(self,problem):
|
||||
# Project ePrimary to bPrimary
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
# Utils used for the data,
|
||||
import numpy as np, matplotlib.pyplot as plt, sys
|
||||
import SimPEG as simpeg
|
||||
import simpegMT as simpegmt
|
||||
|
||||
def getAppRes(MTdata):
|
||||
# Make impedance
|
||||
@@ -78,7 +79,7 @@ def plotMT1DModelData(problem,models,symList=None):
|
||||
else:
|
||||
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
|
||||
# Plot the data and the model
|
||||
colRat = nr/((len(modelList)-2)*1.)
|
||||
colRat = nr/((len(modelList)-1.999)*1.)
|
||||
if colRat > 1.:
|
||||
col = 'k'
|
||||
else:
|
||||
@@ -125,4 +126,43 @@ def plotMT1DModelData(problem,models,symList=None):
|
||||
|
||||
def printTime():
|
||||
import time
|
||||
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
|
||||
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
|
||||
|
||||
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
|
||||
# Find the unique locations
|
||||
# Need to find the locations
|
||||
recData = MTdata.toRecArray()
|
||||
uniLocs = rec2ndarr(np.unique(recData[['x','y','z']])).data
|
||||
mtData1DList = []
|
||||
if 'zxy' in rxType3D:
|
||||
corr = -1 # Shift the data to comply with the quadtrature of the 1d problem
|
||||
else:
|
||||
corr = 1
|
||||
for loc in uniLocs:
|
||||
# Make the receiver list
|
||||
rx1DList = []
|
||||
for rxType in ['z1dr','z1di']:
|
||||
rx1DList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(loc,2).T,rxType))
|
||||
# Source list
|
||||
locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
|
||||
dat1DList = []
|
||||
src1DList = []
|
||||
for src in MTdata.survey.srcList:
|
||||
src1DList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rx1DList,src.freq))
|
||||
for comp in ['r','i']:
|
||||
dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== src.freq].data )
|
||||
|
||||
# Make the survey
|
||||
sur1D = simpegmt.SurveyMT.SurveyMT(src1DList)
|
||||
|
||||
# Make the data
|
||||
dataVec = np.hstack(dat1DList)
|
||||
dat1D = simpegmt.DataMT.DataMT(sur1D,dataVec)
|
||||
sur1D.dobs = dataVec
|
||||
# Need to take MTdata.survey.std and split it as well.
|
||||
std=0.05
|
||||
sur1D.std = np.abs(sur1D.dobs*std) + 0.01*np.linalg.norm(sur1D.dobs)
|
||||
mtData1DList.append(dat1D)
|
||||
|
||||
# Return the the list of data.
|
||||
return mtData1DList
|
||||
@@ -131,9 +131,9 @@ class EDIimporter:
|
||||
|
||||
# Hidden functions
|
||||
def _findLatLong(fileLines):
|
||||
latDMS = np.array(fileLines[_findLine(' LAT=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
|
||||
longDMS = np.array(fileLines[_findLine(' LONG=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
|
||||
elevM = np.array([fileLines[_findLine(' ELEV=',fileLines)[0]].split('=')[1].split()[0]],float)
|
||||
latDMS = np.array(fileLines[_findLine('LAT=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
|
||||
longDMS = np.array(fileLines[_findLine('LONG=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
|
||||
elevM = np.array([fileLines[_findLine('ELEV=',fileLines)[0]].split('=')[1].split()[0]],float)
|
||||
# Convert to D.ddddd values
|
||||
latS = np.sign(latDMS[0])
|
||||
longS = np.sign(longDMS[0])
|
||||
|
||||
Reference in New Issue
Block a user