'k = -1' for ert.simulate

[Coastal0]

When trying to simulate ERT data, I get an odd result.

data = ERT.simulate(mesh = meshERT, res = resBulk, scheme = ertS cheme, verbose=1, returnArray=False) returns the following data;

# a b m n k valid 
1   2   3   4   -1.00000000000000e+000  1
1   2   4   5   -1.00000000000000e+000  1
1   2   5   6   -1.00000000000000e+000  1
...
42  43  45  46  -1.00000000000000e+000  1
43  44  45  46  -1.00000000000000e+000  1

Why would all of my k values be -1?

pyGimli.zip

[Coastal0]

Including noise appears to create an rho_app column... data = ERT.simulate(mesh = meshERT, res=resBulk, scheme=ertScheme, noiseLevel=2, noiseAbs=1e-5)

I guess that's solved....!!

[halbmy]

Unfortunately I cannot follow completely as I do not have the whole script (not included in the zip file) but adding noise to the result should not solve the general problem of synthetic modelling.

The provided scheme file ertScheme already includes a geometric factor of constantly -1, so there is no wonder you also get it in the resulting data. One must think you used k=-1 it for some reason (e.g. inverting resistances but still keeping them positive to use logarithms).

Is the problem (without adding noise) still apparent if you save the schemefile in the generation script with scheme.save('ertScheme', 'a b m n') only?

[halbmy]

mesh = pg.Mesh('meshERT.bms')
print(mesh)
resBulk = np.load('resBulk.npy')
print(resBulk.shape)

yields

Mesh: Nodes: 10627 Cells: 21073 Boundaries: 31699
(78968,)

so the files are not matching.

[Coastal0]

Thanks Thomas. Sorry for the delay, I've finally had a chance to get back to this one. Hopefully I've resolved the mismatch between file and data now, but I still have a k=-1 problem.

mesh+data.zip

print(meshERT)
Mesh: Nodes: 39813 Cells: 78968 Boundaries: 118780

print(resBulk.shape)
(78968,)

If I understand it correctly, all of the cells in meshERThave an associated data from resBulk, so that particular issue should be solved...

Is the problem (without adding noise) still apparent if you save the schemefile in the generation script with scheme.save('ertScheme', 'a b m n') only?

I'm not quite sure I understand the question here, but how can use the ertScheme with a b m n only?

I've used ertScheme = ert.createERTData(sensors, schemeName='dd') to create the array. sensors.zip

It's just occurred to me that the datafile lists the coordinates in x y z despite being a 2D (x z) problem. Could that be related to the geometric factor oddities? data_model.zip Edit; I highly doubt this is the problem, as Florian's flow-modelling example exhibits the same apparent problems. I'm now at a complete and utter loss as to the reason behind the k=-1 issue. :(

[Coastal0]

Just incase it's something simple I've missed, I've copied the major chunks of code below;

Sensor Positions and Creation

sensor_firstX = 50;
sensor_lastX = 250;
sensor_dx= 5;
topo = findTopo(hull) # This is a function which finds the topography at a point (x) along the profile.
sensor_x = np.arange(sensor_firstX,sensor_lastX+sensor_dx,sensor_dx)
sensor_z = np.interp(sensor_x,topo[:,0],topo[:,1])
sensors = np.stack([sensor_x,np.around(sensor_z,2)]).T

ertScheme = ert.createERTData(sensors, schemeName='dd')
ertScheme.save('ertScheme')
print(ertScheme)

ertScheme.txt

Create modelling mesh

meshERT = pg.meshtools.createParaMesh(sensors = ertScheme, quality = 32, paraMaxCellSize=2, paraDepth = 80, paraDX = 0.33)
print('ERT Mesh = ',meshERT)

meshERT.zip

Fill modelling mesh

# Check for negative numbers, set to arbitrary value.
# (This will generally only happen for a bad feflow mesh.)
if any(dInterp <= 1):
    dInterp[dInterp <= 1] = 100

# Convert mass concnetration to water conductivity 
print(dInterp)
k = 0.612 # Linear conversion factor from TDS to EC
sigmaFluid = dInterp/(k*10000) #dInterp (mg/L) to fluid conductivity (S/m)
print(sigmaFluid)
rFluid = 1/sigmaFluid
print(rFluid)

resBulk = petro.resistivityArchie(rFluid, porosity=0.3, m=2, mesh = bPolyMesh, meshI = meshERT, fill=1) 
print(resBulk.shape)

# apply background resistivity model
rho0 = np.zeros(meshERT.cellCount()) + 1000. #Set background resistivity to a number
for c in meshERT.cells():
    if c.center()[1] < 0:
        rho0[c.id()] = 100.
    elif c.center()[1] > 0:
        rho0[c.id()] = 2000.
#pg.show(meshERT,resBulk,colorBar=True,cmap = 'jet',showMesh = True)    

# Apply non-aquifer resistivity
for c in meshERT.cells():
    if c.center()[1] > 0:
        resBulk[c.id()] = 1000.
    elif c.center()[1] < -30:
        resBulk[c.id()] = 20.       

resBulk.zip

Forward model and Invert

data = ert.ERTManager.simulate(mesh = meshERT, res = resBulk, scheme = ertScheme, verbose=1, returnArray=False)
data.save('data')

data.txt

Incase you get a strange desire to see the rest of the code, I've attached it (and the requisite data) below.

data.zip

Workbook_Streamlined_v2.py.txt

[halbmy]

Well I never really worked with the lightweight ERTManager inside pygimli.physics.ert but with the ERTManager from pybert (that should be moved to the first one once it covers everything). By switching to this (importing ERTManager from pybert, and also createData) everything seems to work until the point of computation. I am attaching the script... Workbook_Streamlined_v3TG.py.txt

[Coastal0]

That works perfectly. Thankyou!

So the solution was to use: ertScheme = pb.createData(sensors, schemeName='dd') instead of: ertScheme = ert.createERTData(sensors, schemeName='dd') and ertScheme.set('k', pb.geometricFactor(ertScheme)) ?