Skip to content

Commit

Permalink
Merge branch 'dev' of https://github.com/gimli-org/gimli into dev
Browse files Browse the repository at this point in the history
  • Loading branch information
Thomas Günther committed Jun 19, 2023
2 parents d37a0eb + 94e6e16 commit b6eef3e
Show file tree
Hide file tree
Showing 6 changed files with 239 additions and 226 deletions.
5 changes: 3 additions & 2 deletions CMakeLists.txt
View file
Expand Up @@ -67,8 +67,9 @@ if (CMAKE_COMPILER_IS_GNUCXX OR CMAKE_COMPILER_IS_CLANGXX OR ${CMAKE_CXX_COMPILE
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-int-in-bool-context")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unused-local-typedefs") # pg build
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated-declarations") # loot of boost complaints about bind
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unknown-warning-option") # on appleclang
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-undefined-var-template") # on appleclang
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unknown-warning-option") # on apple clang
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-undefined-var-template") # on apple clang
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-overloaded-virtual")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DBOOST_BIND_GLOBAL_PLACEHOLDERS")

set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -ansi -pedantic -fno-omit-frame-pointer -Wall")
Expand Down
11 changes: 8 additions & 3 deletions doc/tutorials/3_inversion/plot_8-regionWise.py
View file
Expand Up @@ -108,7 +108,7 @@

mgr = ert.ERTManager(data, verbose=True)
mgr.setMesh(mesh) # use this mesh for all subsequent runs
mgr.invert()
mgr.invert(maxIter=1)
# mgr.invert(mesh=mesh) would only temporally use the mesh

# %%%
Expand Down Expand Up @@ -136,13 +136,18 @@
# %%%
# Apparently, the two regions are already decoupled from each other which
# makes sense. Let us look in detail at the water cells by extracting the
# water body.
# water body.
# Note. The manager class performs a model value permutation to fit
# the parametric mesh cell. So if you want to relate model values to the input
# mesh, you need to use the unpermutated model values directly from the
# inversion framework instance: `mgr.fw.model``
#

water = mesh.createSubMesh(mesh.cells(mesh.cellMarkers() == 3))
resWater = mgr.model[len(mgr.model)-water.cellCount():]
resWater = mgr.fw.model[len(mgr.model)-water.cellCount():]
ax, cb = pg.show(water, resWater)


# %%%
# Apparently, all values are below the expected 22.5\ $\Omega$\ m
# and some are implausibly low. Therefore we should try to limit them.
Expand Down
5 changes: 5 additions & 0 deletions pygimli/frameworks/methodManager.py
View file
Expand Up @@ -709,6 +709,11 @@ def __init__(self, **kwargs):
super(MeshMethodManager, self).__init__(**kwargs)
self.mesh = None

@property
def model(self):
"""Inversion model."""
return self.paraModel(self.fw.model)

@property
def paraDomain(self):
"""Parameter (inversion) domain mesh."""
Expand Down
4 changes: 2 additions & 2 deletions pygimli/frameworks/modelling.py
View file
Expand Up @@ -740,10 +740,10 @@ def drawModel(self, ax, model, **kwargs):
# TODO needs to be checked if mapping is always ok (region example)
# is (len(model) == self.paraDomain.cellCount() or \
if hasattr(model, "isParaModel") and model.isParaModel is False:
pg._y(model.isParaModel)
# pg._y(model.isParaModel)
mod = self.paraModel(model)
elif hasattr(model, "isParaModel") and model.isParaModel is True:
pg._g(model.isParaModel)
# pg._g(model.isParaModel)
mod = model
elif len(model) == self.paraDomain.nodeCount():
# why nodeCount? a field as model result?
Expand Down
120 changes: 60 additions & 60 deletions pygimli/physics/em/hemmodelling.py
View file
Expand Up @@ -72,7 +72,7 @@ def __init__(self, nlay, height, f=None, r=None, **kwargs):
if self.f is None:
self.f = self.fdefault
if isinstance(r, float) or isinstance(r, int):
self.r = np.ones_like(f, dtype=np.float) * r
self.r = np.ones_like(f, dtype=float) * r
else:
if len(r) == len(self.f):
self.r = r
Expand All @@ -82,7 +82,7 @@ def __init__(self, nlay, height, f=None, r=None, **kwargs):
self.r = self.rdefault

self.wem = (2.0 * pi * self.f) ** 2 * self.ep0 * self.mu0
self.iwm = np.complex(0, 1) * 2.0 * pi * self.f * self.mu0
self.iwm = 1.0j * 2.0 * pi * self.f * self.mu0
mesh = pg.meshtools.createMesh1DBlock(nlay)
super().__init__()
self.setMesh(mesh)
Expand All @@ -100,40 +100,40 @@ def response_mt(self, par, i=0):

def calc_forward(self, x, h, rho, d, epr, mur, quasistatic=False):
"""Calculate forward response."""
field = np.zeros((self.f.size, x.size), np.complex)
field = np.zeros((self.f.size, x.size), complex)
# Forward calculation for background model
if d.size:
for m in range(x.size):
field[:, m] = self.vmd_hem(np.array([h[m]], np.float),
field[:, m] = self.vmd_hem(np.array([h[m]], float),
rho[:, m], d[:, m], epr[:, m],
mur[:, m], quasistatic).T[:, 0]
else:
for n in range(self.f.size):
for m in range(x.size):
field[n, m] = self.vmd_hem(np.array([h[n, m]], np.float),
np.array([rho[n, m]], np.float),
field[n, m] = self.vmd_hem(np.array([h[n, m]], float),
np.array([rho[n, m]], float),
d,
np.array([epr[n, m]], np.float),
np.array([mur[n, m]], np.float),
np.array([epr[n, m]], float),
np.array([mur[n, m]], float),
quasistatic).T[n, 0]
return field

def downward(self, rho, d, z, epr, mur, lam):
"""Downward continuation of fields."""
nl = rho.size
alpha = np.zeros((nl, lam.shape[1], self.f.size), np.complex)
b = np.zeros((nl, lam.shape[1], self.f.size), np.complex)
aa = np.zeros((nl, lam.shape[1], self.f.size), np.complex)
aap = np.zeros((nl, lam.shape[1], self.f.size), np.complex)
alpha = np.zeros((nl, lam.shape[1], self.f.size), complex)
b = np.zeros((nl, lam.shape[1], self.f.size), complex)
aa = np.zeros((nl, lam.shape[1], self.f.size), complex)
aap = np.zeros((nl, lam.shape[1], self.f.size), complex)
rho = rho[:, np.newaxis, np.newaxis] * np.ones(
(rho.size, lam.shape[1], self.f.size), np.float)
(rho.size, lam.shape[1], self.f.size), float)
d = d[:, np.newaxis, np.newaxis] * np.ones(
(d.size, lam.shape[1], self.f.size), np.float)
(d.size, lam.shape[1], self.f.size), float)
h = np.insert(np.cumsum(d[:, 0, 0]), 0, 0)
epr = epr[:, np.newaxis, np.newaxis] * np.ones(
(epr.size, lam.shape[1], self.f.size), np.float)
(epr.size, lam.shape[1], self.f.size), float)
mur = mur[:, np.newaxis, np.newaxis] * np.ones(
(mur.size, lam.shape[1], self.f.size), np.float)
(mur.size, lam.shape[1], self.f.size), float)
lam = np.tile(lam, (nl, 1, 1))
# progression constant
alpha = np.sqrt(lam ** 2 - np.tile(self.wem, (nl, lam.shape[1], 1)) *
Expand Down Expand Up @@ -211,24 +211,24 @@ def vmd_hem(self, h, rho, d, epr=1., mur=1., quasistatic=False):
"""
# filter coefficients
if isinstance(epr, float):
epr = np.ones((len(rho),), np.float)*epr
epr = np.ones((len(rho),), float)*epr
if isinstance(mur, float):
mur = np.ones((len(rho),), np.float)*mur
mur = np.ones((len(rho),), float)*mur
fc0, nc, nc0 = hankelfc(3)
fc1, nc, nc0 = hankelfc(4)
# allocate arrays
nf = len(self.f)
lam = np.zeros((1, nc, nf), np.float)
alpha0 = np.zeros((1, nc, nf), np.complex)
delta0 = np.zeros((1, nc, nf), np.complex)
delta1 = np.zeros((1, nc, nf), np.complex)
delta2 = np.zeros((1, nc, nf), np.complex)
delta3 = np.zeros((1, nc, nf), np.complex)
aux0 = np.zeros((1, nf), np.complex)
aux1 = np.zeros((1, nf), np.complex)
aux2 = np.zeros((1, nf), np.complex)
aux3 = np.zeros((1, nf), np.complex)
Z = np.zeros(nf, np.complex)
lam = np.zeros((1, nc, nf), float)
alpha0 = np.zeros((1, nc, nf), complex)
delta0 = np.zeros((1, nc, nf), complex)
delta1 = np.zeros((1, nc, nf), complex)
delta2 = np.zeros((1, nc, nf), complex)
delta3 = np.zeros((1, nc, nf), complex)
aux0 = np.zeros((1, nf), complex)
aux1 = np.zeros((1, nf), complex)
aux2 = np.zeros((1, nf), complex)
aux3 = np.zeros((1, nf), complex)
Z = np.zeros(nf, complex)
# r0
r0 = np.copy(self.r)
# determine optimum r0 (shift nodes) for f > 1e4 and h > 100
Expand All @@ -243,13 +243,13 @@ def vmd_hem(self, h, rho, d, epr=1., mur=1., quasistatic=False):
r0[index] = self.c0 / (2.0 * np.pi * self.f[index]) * 10.0 ** (
(opt + 0.5 - nc0) / 10.0)
# Wave numbers
n = np.arange(nc0 - nc, nc0, 1, np.float)
n = np.arange(nc0 - nc, nc0, 1, float)
q = 0.1 * np.log(10)
lam = np.reshape(np.exp(-n[np.newaxis, :, np.newaxis] * q) /
r0[np.newaxis, np.newaxis, :], (-1, nc, r0.size))
# (1, 100, nfreq)
# wave number in air, quasistationary approximation
alpha0 = np.copy(lam) * np.complex(1, 0) # (1, 100, nfreq)
alpha0 = np.copy(lam) * complex(1, 0) # (1, 100, nfreq)
# wave number in air, full solution for f > 1e4
if quasistatic:
index = np.zeros(self.f.shape, np.bool)
Expand All @@ -271,18 +271,18 @@ def vmd_hem(self, h, rho, d, epr=1., mur=1., quasistatic=False):
# quasistationary approximation
aux0 = np.sum(delta0 * lam ** 3 / alpha0 *
np.tile(fc0[::-1].T[:, :, np.newaxis],
(1, 1, self.f.size)), 1, np.complex) / r0
(1, 1, self.f.size)), 1, complex) / r0
# full solution, partial integration
if np.any(index):
aux1 = np.sum(delta1 * lam ** 3 *
np.tile(fc0[::-1].T[:, :, np.newaxis],
(1, 1, self.f.size)), 1, np.complex) / r0
(1, 1, self.f.size)), 1, complex) / r0
aux2 = np.sum(
delta2 * lam * np.tile(fc0[::-1].T[:, :, np.newaxis],
(1, 1, self.f.size)), 1, np.complex)/r0
(1, 1, self.f.size)), 1, complex)/r0
aux3 = np.sum(delta3 * lam ** 2 *
np.tile(fc1[::-1].T[:, :, np.newaxis],
(1, 1, self.f.size)), 1, np.complex) / r0
(1, 1, self.f.size)), 1, complex) / r0
# normed secondary field
# quasistationary approximation
Z = self.r ** 3 * aux0 * self.scaling
Expand All @@ -298,14 +298,14 @@ def vmd_total_Ef(self, h, z, rho, d, epr, mur, tm):
# only halfspace
# Filter coefficients
fc1, nc, nc0 = hankelfc(4)
lam = np.zeros((1, nc, self.f.size), np.float)
alpha0 = np.zeros((1, nc, self.f.size), np.complex)
delta = np.zeros((1, nc, self.f.size), np.complex)
aux = np.zeros((1, self.f.size), np.complex)
Ef = np.zeros(self.f.size, np.complex)
lam = np.zeros((1, nc, self.f.size), float)
alpha0 = np.zeros((1, nc, self.f.size), complex)
delta = np.zeros((1, nc, self.f.size), complex)
aux = np.zeros((1, self.f.size), complex)
Ef = np.zeros(self.f.size, complex)
r0 = np.copy(self.r)
# wave numbers
n = np.arange(nc0 - nc, nc0, 1, np.float)
n = np.arange(nc0 - nc, nc0, 1, float)
q = 0.1 * np.log(10)
lam = np.reshape(np.exp(-n[np.newaxis, :, np.newaxis] * q) /
r0[np.newaxis, np.newaxis, :], (-1, nc, r0.size))
Expand All @@ -321,7 +321,7 @@ def vmd_total_Ef(self, h, z, rho, d, epr, mur, tm):
# quasistationary approximation
aux = np.sum(delta*lam**2*aa*np.tile(fc1[::-1].T[:, :, np.newaxis],
(1, 1, self.f.size)), 1,
np.complex) / r0 # (1, nfreq)
complex) / r0 # (1, nfreq)
# absolute fields, full solution
Ef = -tm * self.iwm / (4.0 * np.pi) * aux
return Ef
Expand Down Expand Up @@ -351,9 +351,9 @@ def hankelfc(order):
-5.35535069e-5, 3.37899801e-5, -2.13200365e-5, 1.34520337e-5,
-8.48765949e-6, 5.35535110e-6, -3.37899811e-6, 2.13200368e-6,
-1.34520338e-6, 8.48765951e-7, -5.35535110e-7, 3.37899811e-7],
np.float)
nc = np.int(80)
nc0 = np.int(40)
float)
nc = int(80)
nc0 = int(40)
elif order == 2: # cos
fc = np.array([
1.63740363e-7, 1.83719709e-7, 2.06136904e-7, 2.31289411e-7,
Expand Down Expand Up @@ -397,9 +397,9 @@ def hankelfc(order):
2.73337984e-5, -1.72464607e-5, 1.08817810e-5, -6.86593962e-6,
4.33211503e-6, -2.73337979e-6, 1.72464606e-6, -1.08817810e-6,
6.86593961e-7, -4.33211503e-7, 2.73337979e-7, -1.72464606e-7],
np.float)
nc = np.int(164)
nc0 = np.int(122)
float)
nc = int(164)
nc0 = int(122)
elif order == 3: # J0
fc = np.array([
2.89878288e-7, 3.64935144e-7, 4.59426126e-7, 5.78383226e-7,
Expand Down Expand Up @@ -427,9 +427,9 @@ def hankelfc(order):
9.89181741e-6, -6.24131160e-6, 3.93800058e-6, -2.48471018e-6,
1.56774609e-6, -9.89180896e-7, 6.24130948e-7, -3.93800005e-7,
2.48471005e-7, -1.56774605e-7, 9.89180888e-8, -6.24130946e-8],
np.float)
nc = np.int(100)
nc0 = np.int(60)
float)
nc = int(100)
nc0 = int(60)
elif order == 4: # J1
fc = np.array([
1.84909557e-13, 2.85321327e-13, 4.64471808e-13, 7.16694771e-13,
Expand Down Expand Up @@ -457,9 +457,9 @@ def hankelfc(order):
-8.64396364e-5, 5.45397224e-5, -3.44122382e-5, 2.17126544e-5,
-1.36997587e-5, 8.64396338e-6, -5.45397218e-6, 3.44122380e-6,
-2.17126543e-6, 1.36997587e-6, -8.64396337e-7, 5.45397218e-7],
np.float)
nc = np.int(100)
nc0 = np.int(60)
float)
nc = int(100)
nc0 = int(60)
return (np.reshape(fc, (-1, 1)), nc, nc0) # (100,) -> (100, 1)


Expand Down Expand Up @@ -494,9 +494,9 @@ def __init__(self, nlay, height=1, **kwargs):

def response(self, par):
"""Response vector as combined in-phase and out-phase data."""
thk = np.asarray(par[:self.nlay-1], dtype=np.float)
res = np.asarray(par[self.nlay-1:2*self.nlay-1], dtype=np.float)
mur = np.asarray(par[2*self.nlay-1:3*self.nlay-1], dtype=np.float) + 1
thk = np.asarray(par[:self.nlay-1], dtype=float)
res = np.asarray(par[self.nlay-1:2*self.nlay-1], dtype=float)
mur = np.asarray(par[2*self.nlay-1:3*self.nlay-1], dtype=float) + 1
ip, op = self.vmd_hem(self.height, rho=res, d=thk, mur=mur)
return pg.cat(ip, op)

Expand Down Expand Up @@ -605,9 +605,9 @@ def createJacobian(self, model):

if __name__ == '__main__':
numlay = 3
height = np.float(30.0)
resistivity = np.array([1000.0, 100.0, 1000.0], np.float)
thickness = np.array([22.0, 29.0], np.float)
height = float(30.0)
resistivity = np.array([1000.0, 100.0, 1000.0], float)
thickness = np.array([22.0, 29.0], float)
f = HEMmodelling(numlay, height, r=10) # frequency, separation)
IP, OP = f.vmd_hem(height, resistivity, thickness)
print(IP, OP)

0 comments on commit b6eef3e

Please sign in to comment.