Open
Description
Problem description
I am trying to make the LCI code of TDEM run properly but the results are not good I tried everything but am not getting the best results. Please help me regarding this.
Please describe your issue here.
Your environment
Please provide the output of print(pygimli.Report()) here. If that does not
work, please give provide some additional information on your:
Date: Wed Apr 03 22:46:02 2024 China Standard Time
OS : Windows
CPU(s) : 32
Machine : AMD64
Architecture : 64bit
RAM : 31.7 GiB
Environment : Jupyter
Python 3.9.18 | packaged by conda-forge | (main, Dec 23 2023, 16:29:04) [MSC
v.1929 64 bit (AMD64)]
pygimli : 1.5.0.post1
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.8.0
scipy : 1.12.0
IPython : 8.15.0
meshio : 5.3.5
tetgen : 0.6.4
pyvista : 0.43.4
Intel(R) oneAPI Math Kernel Library Version 2023.2-Product Build 20230613
for Intel(R) 64 architecture applications
import matplotlib.pyplot as plt
import numpy as np
import pygimli as pg
from pygimli.viewer.mpl import drawModel1D, showStitchedModels
from pygimli.physics.em import VMDTimeDomainModelling
class TDEM2dFOP(pg.core.ModellingBase):
"""TDEM 2d-LCI modelling class based on BlockMatrices."""
def __init__(self, x, times, txArea, nlay=2, verbose=False):
super(TDEM2dFOP, self).__init__(verbose)
self.nlay = nlay
self.header = {}
self.pos = None
self.z = None
self.topo = None
self.nx = len(x)
self.nt = len(times)
npar = 2 * nlay - 1 # number of parameters per sounding
self.mesh1d = pg.meshtools.createMesh1D(self.nx, npar)
self.mesh_ = pg.meshtools.createMesh1D(self.nx, npar)
self.setMesh(self.mesh_)
self.fop = pg.physics.em.VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=txArea)
self.jacobian = pg.matrix.BlockMatrix()
self.fop_list = []
for _ in range(self.nx):
self.fop_list.append(pg.physics.em.VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=txArea))
n = self.jacobian.addMatrix(self.fop_list[-1].jacobian())
self.jacobian.addMatrixEntry(n, self.nt * _, npar * _)
self.jacobian.recalcMatrixSize()
print(self.jacobian.rows(), self.jacobian.cols())
self.setJacobian(self.jacobian)
def response(self, model):
"""Combine forward responses of all soundings."""
model_array = np.asarray(model).reshape((self.nx, self.nlay * 2 - 1))
response = pg.Vector(0)
for model_i in model_array:
response = pg.cat(response, self.fop.response(model_i))
return response
def createJacobian(self, model):
"""Create Jacobian matrix by creating individual Jacobians."""
model_array = np.asarray(model).reshape((self.nx, self.nlay * 2 - 1))
for i, model_i in enumerate(model_array):
self.fop_list[i].createJacobian(model_i)
if __name__ == "__main__":
# Define the model parameters
num_positions = 20 # number of sampling positions
layer_resistivities = [1000, 100, 500, 20]
num_layers = len(layer_resistivities) # number of layers
resistivities = np.ones((num_positions, num_layers)) * layer_resistivities
x = np.linspace(0, 20, num_positions)[:, None]
nx = len(x)
thickness_1 = 5 + 0.2 * np.sin(x * np.pi * 2)
thickness_2 = 15 + np.sin(x * np.pi * 2)
thickness_3 = 10 + 0.3 * np.sin(x * np.pi * 2)
# Time-domain parameters
num_time_steps = 64
times = np.logspace(np.log10(1e-4), np.log10(1e-1), num_time_steps)
transmitter_area = 62500
true_model = []
for i in range(num_positions):
model_i = np.array([thickness_1[i][0], thickness_2[i][0], thickness_3[i][0]] + layer_resistivities) # True model
print(model_i)
true_model.append(model_i)
try:
tdem2d_fop = TDEM2dFOP(x, times, transmitter_area, nlay=num_layers)
data_tem = tdem2d_fop.response(true_model)
error_tem_abs = 1. # absolute error (ppm of primary signal)
error_tem = np.ones_like(data_tem) * error_tem_abs
# 2D LCI inversion
# Transformations
trans_data = pg.trans.Trans()
trans_thickness = pg.trans.TransLogLU(0.1, 50.)
trans_resistivity = pg.trans.TransLogLU(10., 2000.)
for i in range(num_layers - 1):
tdem2d_fop.region(i).setTransModel(trans_thickness)
for i in range(num_layers - 1, num_layers * 2 - 1):
tdem2d_fop.region(i).setTransModel(trans_resistivity)
# Set constraints
tdem2d_fop.region(0).setConstraintType(10)
tdem2d_fop.region(1).setConstraintType(10)
# Generate starting model by repetition
starting_model = pg.Vector([10] * (num_layers - 1) + [100] * num_layers)
model = np.repeat(starting_model, len(x))
inversion = pg.core.Inversion(data_tem, tdem2d_fop, trans_data)
inversion.setAbsoluteError(error_tem)
inversion.setLambda(500)
inversion.setModel(model)
inversion.setReferenceModel(model)
estimated_model = inversion.run()
estimated_model_array = np.asarray(estimated_model).reshape((nx, num_layers * 2 - 1))
estimated_model_list = [] # empty array to store the estimated model
for model_i in estimated_model_array:
print(model_i)
estimated_model_list.append(model_i)
# Plotting the results
fig, ax = pg.plt.subplots(figsize=(12, 6), ncols=2)
showStitchedModels(true_model, ax=ax[0], x=None, cMap='Spectral_r', logScale=True, title='True model (Ohm.m)')
ax[0].set_ylim(-50, 0)
ax[0].set_xlabel('Distance (m)')
ax[0].set_ylabel('Depth (m)')
ax[0].set_title('True model')
ax[0].legend(loc='lower right', ncol=2, fontsize=10, frameon=False)
ax[0].grid(True, alpha=0.5)
ax[0].set_axisbelow(True)
showStitchedModels(estimated_model_list, ax=ax[1], x=None, cMap='Spectral_r', logScale=True, title='Estimated model (Ohm.m)')
ax[1].set_ylim(-50, 0)
ax[1].set_xlabel('Distance (m)')
ax[1].set_title('Estimated model')
ax[1].legend(loc='lower right', ncol=2, fontsize=10, frameon=False)
ax[1].grid(True, alpha=0.5)
ax[1].set_axisbelow(True)
plt.show()
except Exception as e:
print(f"An error occurred: {e}")
Metadata
Metadata
Assignees
Labels
No labels