import pygimli as pg
import matplotlib.pyplot as plt
from pybert.manager import Resistivity
from pygimli.utils.base import rmswitherr
# RhoAC = RhoDC * (1-m) ==> dRhoa / m = dRhoa/dRho * dRho/m = JDC * (-RhoDC)
[docs]class DCIPMModelling(pg.ModellingBase):
""" DC/IP modelling class using an (FD-based) approach """
def __init__(self, f, mesh, rho, verbose=False):
""" init class with DC forward operator and resistivity vector """
super().__init__(verbose=verbose)
self.setMesh(mesh)
self.f = f
self.rhoDC = rho # DC resistivity
self.drhodm = -self.rhoDC
self.rhoaDC = f.response(self.rhoDC)
self.dmdrhoa = -1.0 / self.rhoaDC
self.J = pg.matrix.MultLeftRightMatrix(f.jacobian(), self.dmdrhoa,
self.drhodm)
self.setJacobian(self.J)
self.fullJacocian = False
[docs] def response(self, m):
""" return forward response as function of chargeability model """
self.rhoaAC = self.f.response(self.rhoDC * (1. - m))
if self.fullJacocian:
self.dmdrhoa = -1.0 / self.rhoaAC
return pg.abs(1.0 - self.rhoaAC / self.rhoaDC)
[docs] def createJacobian(self, model):
""" create jacobian matrix using unchanged DC jacobian and m model """
if self.fullJacocian:
self.rhoAC = self.rhoDC * (1 - model)
self.J.right = -self.rhoAC
self.f.createJacobian(self.rhoAC)
if __name__ == "__main__":
zw = 0.5
dc = Resistivity('Onneslov.data')
print("DC min/max = ", min(dc.data('ip')), max(dc.data('ip')))
dc.setMeshPot()
dc.fop.setThreadCount(8)
if 0:
dc.invert(zweight=zw) # , recalcJacobian=False, maxIter=2)
dc.showResultAndFit(cMin=0.5, cMax=500)
dc.figs['resultFit'].savefig('resultFit.pdf')
dc.mesh.save('dcmesh.bms')
dc.resistivity.save('res.vec')
else:
dc.resistivity = pg.RVector('resistivity.vector')
dc.response = dc.fop(dc.resistivity)
dc.fop.createJacobian(dc.resistivity)
dc.paraDomain = dc.fop.regionManager().paraDomain()
# %%
# iperr = np.array(dc.error)
ma = dc.data('ip') * 1e-3 # mrad to rad
iperr = pg.RVector(dc.data.size(), 0.01)
if 1:
mmin, mmax = 1e-3, 0.2
print('discarding min/max', sum(ma < mmin), sum(ma > mmax))
iperr[ma < mmin] = 1e5
ma[ma < mmin] = mmin
iperr[ma > mmax] = 1e5
ma[ma > mmax] = mmax
# %% set up forward operator
fIP = DCIPMModelling(dc.fop, dc.mesh, dc.resistivity)
fIP.region(1).setBackground(True)
fIP.region(2).setConstraintType(1)
fIP.region(2).setZWeight(zw)
fIP.createRefinedForwardMesh(True)
tD, tM = pg.RTransLog(), pg.RTransLogLU(0, 0.95)
INV = pg.RInversion(ma, fIP, tD, tM, True, False)
mstart = pg.RVector(len(dc.resistivity), pg.median(ma)) # 10 mV/V
INV.setModel(mstart)
INV.setAbsoluteError(iperr)
INV.setLambda(10)
INV.setRobustData(True)
dc.m = INV.run()
dc.m.save("m.vec")
dc.paraDomain = fIP.regionManager().paraDomain()
cMin, cMax = 10, 300 # max(ma)*1e3
# %%
respma = INV.response() # AC/DC dc.data('rhoa')
print(max(ma), max(respma), max(dc.m))
myrms = rmswitherr(ma, respma, INV.error(), 0.4)
print("RMS=", myrms)
fig, ax = plt.subplots(nrows=3, figsize=(8, 12))
# dc.showModel(ax=ax[2], vals=dc.m*1e3, cMin=cMin, cMax=cMax, logScale=False)
pg.show(dc.paraDomain, dc.m*1e3, ax=ax[2], cMin=cMin, cMax=cMax,
logScale=True, colorBar=True)
for i, mdata in enumerate([ma, respma]):
dc.showData(ax=ax[i], vals=mdata*1e3, cMin=cMin, cMax=cMax, logScale=False,
colorBar=False)
# %%
fig.savefig('outMaM.pdf', bbox_inches='tight')
