From 291cce9d5227b2e83564ea44dc80b8c729aa7fe0 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 13 Nov 2015 12:49:19 -0800 Subject: [PATCH] Finish forward modeling in integral Start test function. --- simpegPF/Dev/MAG3D_fwr_Driver.py | 5 +++-- simpegPF/Dev/MAG3D_fwr_Test.py | 16 +++++++--------- simpegPF/Dev/fwr_MAG_obs.py | 7 ++++++- simpegPF/Magnetics.py | 4 ++++ 4 files changed, 20 insertions(+), 12 deletions(-) diff --git a/simpegPF/Dev/MAG3D_fwr_Driver.py b/simpegPF/Dev/MAG3D_fwr_Driver.py index 7e1ff802..f982563c 100644 --- a/simpegPF/Dev/MAG3D_fwr_Driver.py +++ b/simpegPF/Dev/MAG3D_fwr_Driver.py @@ -1,6 +1,6 @@ import os -home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\Research\SimPeg' +home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev' inpfile = 'MAG3Cfwr.inp' @@ -28,9 +28,10 @@ model = Utils.meshutils.readUBCTensorModel(modfile,mesh) [B,M,dobs] = BaseMag.readUBCmagObs(obsfile) rxLoc = dobs[:,0:3] +ndata = rxLoc.shape[0] # Compute forward model using integral equation -d = fwr_MAG_obs(xn,yn,zn,B,M,rxLoc,model) +d = fwr_MAG_obs(mesh,B,M,rxLoc,model) # Form data object with coordinates and write to file data = np.c_[rxLoc , d , np.zeros((ndata,1))] diff --git a/simpegPF/Dev/MAG3D_fwr_Test.py b/simpegPF/Dev/MAG3D_fwr_Test.py index 03665e86..aa97e332 100644 --- a/simpegPF/Dev/MAG3D_fwr_Test.py +++ b/simpegPF/Dev/MAG3D_fwr_Test.py @@ -1,6 +1,6 @@ import os -home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\Research\SimPeg' +home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev' os.chdir(home_dir) @@ -15,12 +15,10 @@ from fwr_MAG_obs import fwr_MAG_obs #%% Create survey -B[0] = 90. -B[1] = 00. -B[2] = 50000. +B = np.array(([90.,0.,50000.])) -M[0] = 90. -M[1] = 0. +M = np.array(([90.,0.])) + # # Or create juste a plane grid xr = np.linspace(-1./2., 1./2., 10) @@ -55,13 +53,13 @@ model = np.ones(mcell)*chibkg model[sph_ind] = chiblk #%% Forward mode ldata -d = fwr_MAG_obs(xn,yn,zn,B,M,rxLoc,model) +d = fwr_MAG_obs(mesh,B,M,rxLoc,model) #%% Get the analystical answer and compute the residual bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],.25,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary') r_Bz = mkvc(d) - bza -lrl = sum( dBz**2 ) **0.5 +lrl = sum( r_Bz**2 ) **0.5 print "Residual between analytical and integral= " + str(lrl) #%% Plot fields @@ -95,6 +93,6 @@ plt.scatter(X,Y, c='k', s=5) #%% Compare fields plt.subplot(122) -plt.imshow(np.reshape(dBz,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()]) +plt.imshow(np.reshape(r_Bz,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()]) plt.contour(X,Y, np.reshape(dBz,X.shape),10) plt.scatter(X,Y, c='k', s=5) \ No newline at end of file diff --git a/simpegPF/Dev/fwr_MAG_obs.py b/simpegPF/Dev/fwr_MAG_obs.py index 6b65ffe5..623ea89a 100644 --- a/simpegPF/Dev/fwr_MAG_obs.py +++ b/simpegPF/Dev/fwr_MAG_obs.py @@ -1,4 +1,4 @@ -def fwr_MAG_obs(xn,yn,zn,B,M,rxLoc,model): +def fwr_MAG_obs(mesh,B,M,rxLoc,model): """ Forward model magnetic data using integral equation @@ -21,6 +21,11 @@ def fwr_MAG_obs(xn,yn,zn,B,M,rxLoc,model): from SimPEG import np, Utils, sp, mkvc from get_T_mat import get_T_mat + + xn = mesh.vectorNx; + yn = mesh.vectorNy; + zn = mesh.vectorNz; + mcell = (len(xn)-1) * (len(yn)-1) * (len(zn)-1) ndata = rxLoc.shape[0] diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index 9d8589e9..dc96baa8 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -3,6 +3,10 @@ import BaseMag from scipy.constants import mu_0 from MagAnalytics import spheremodel, CongruousMagBC +class MagneticIntegral(Problem.BaseProblem): + """ + approach using IE + """ class MagneticsDiffSecondary(Problem.BaseProblem):