01_Sim.py

#%%
import matplotlib.pyplot as plt
import numpy as np

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
import pandas as pd

world = mt.createWorld(start=[0, 0], end=[71, -10], layers=[-1.5],
                       worldMarker=True)

# block = mt.createCircle(pos=[30, -8.], radius=[4, 1], marker=2,
#                         boundaryMarker=10, area=0.1)  

poly = mt.createPolygon([(7,-3), (13,-3),(38,-6.5), (30,-6.5)
                          ], isClosed=True,
                          addNodes=3, interpolate='spline', marker=2)     

water_level = mt.createPolygon([(0,-8), (30,-7.5), (50,-8), (60,-7.5),(71,-7.8),(71,-10),
                                (0,-10)
                          ], isClosed=True,
                          addNodes=0, interpolate='spline', marker=2)

# geom = world + block + poly + water_level
geom = world + poly +  water_level

pg.show(geom)
# %%

scheme = ert.createERTData(elecs=np.linspace(start=0, stop=71, num=72),
                           schemeName='dd')
for p in scheme.sensors():
    geom.createNode(p)
    geom.createNode(p - [0, 0.1])

# Create a mesh for the finite element modelling with appropriate mesh quality.
mesh = mt.createMesh(geom, quality=34,area =0.5)


mesh.save('./A_Simresults/Simulation_mesh')

# Create a map to set resistivity values in the appropriate regions
# [[regionNumber, resistivity], [regionNumber, resistivity], [...]
rhomap = [[0,35],
          [1, 80.],
          [2, 10.]]

# Take a look at the mesh and the resistivity distribution
pg.show(mesh, data=rhomap, label=pg.unit('res'), showMesh=True)
# %%
# inverse_par = pd.Series({'plotmin':1,'plotmax':100,'zmin':0,'zmax':29.2,'dx':0.25,'lambda_reg':10,\
#     'paraDX':0.3,'paraMaxCellSize':0.5,'paraDepth':30,'quality':34})
    
    
data = ert.simulate(mesh, scheme=scheme, res=rhomap, noiseLevel=1,
                    noiseAbs=1e-6, seed=1337)

plt.close('all')
pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
pg.info('Simulated data', data)
pg.info('The data contains:', data.dataMap().keys())

pg.info('Simulated rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))
pg.info('Selected data noise %(min/max)', min(data['err'])*100, max(data['err'])*100)

data.remove(data['rhoa'] < 0)
pg.info('Filtered rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))

# You can save the data for further use
data.save('./A_Simresults/Simulated_resistance.dat')


# You can take a look at the data
ert.show(data)
pg.show(mesh, rhomap, hold=True, cMap="Spectral_r", logScale=True,
        orientation="vertical", cMin=0, cMax=100)