Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Create a 3D layer model #695

Open
zhichaoluan opened this issue Apr 11, 2024 · 11 comments
Open

Create a 3D layer model #695

zhichaoluan opened this issue Apr 11, 2024 · 11 comments

Comments

@zhichaoluan
Copy link

from pygimli.physics import ert
import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
import pygimli.meshtools as mt

定义模型层和标记

top_layer = mt.createWorld(start=[-10, -10, 0], end=[140, 10, -10], marker=1)
middle_layer = mt.createWorld(start=[-10, -10, -10], end=[140, 10, -20], marker=2)
bottom_layer = mt.createWorld(start=[-10, -10, -20], end=[140, 10, -40], marker=3)
cube_region = mt.createCube(size=[128, 0.4, 40], pos=[128/2, 0, -20], marker=4)

合并模型区域

start_model = mt.merge([top_layer, middle_layer, bottom_layer, cube_region])

展示初始模型

pg.show(start_model, alpha=0.3, markers=True)

加载数据

filename = "cc.shm"
cc.zip

shm = pg.DataContainerERT(filename)
data = ert.load(filename)
data['k'] = ert.geometricFactors(data)

计算视电阻率

data['rhoa'] = data['k'] * data['u'] / data['i'] * 1000
data.remove(data['rhoa'] < 0)
data['err'] = ert.estimateError(data, relativeError=0.02)

在模型中添加传感器位置

for s in shm.sensors():
start_model.createNode(s)
for s in shm.sensorPositions():
start_model.createNode(s - [0, 0, 1e-2 / 2])

创建反演网格

inversion_mesh = mt.createMesh(start_model, quality=1.4)
pg.show(inversion_mesh, markers=True, showMesh=True)

配置ERT管理器

mgr = ert.ERTManager()
mgr.setData(data)
mgr.setMesh(inversion_mesh)

设置正则化参数

mgr.inv.setRegularization(1, fix=100) # 顶层
mgr.inv.setRegularization(2, fix=200) # 中层
mgr.inv.setRegularization(3, fix=300) # 底层
mgr.inv.setRegularization(4, startModel=1e4, limits=[100, 1e4 + 20]) # 立方体区域

执行反演

mgr.invert(lam=10, zWeight=0.1, verbose=True)

保存结果并展示误差

mgr.saveResult()
mgr.showMisfit()

I want to create a layer model and assign different resistivity values to each for constraint inversion, but I have some problems to help me solve.

Traceback (most recent call last):
File "D:/wangkangbo/test/11111.py", line 36, in
mgr.invert(lam=10,zWeight=0.1,verbose=True)
File "D:\Anaconda3\envs\pg1\lib\site-packages\pygimli\frameworks\methodManager.py", line 777, in invert
if self.fop.mesh() is None:
File "D:\Anaconda3\envs\pg1\lib\site-packages\pygimli\frameworks\modelling.py", line 665, in mesh
self._applyRegionProperties()
File "D:\Anaconda3\envs\pg1\lib\site-packages\pygimli\frameworks\modelling.py", line 325, in _applyRegionProperties
if rMgr.region(rID).fixValue() != vals['fix']:
RuntimeError: C:/msys64/home/halbm/gimli/gimli/core/src/regionManager.cpp:593 GIMLI::Region* GIMLI::RegionManager::region(GIMLI::SIndex) no region with marker 1

@halbmy
Copy link
Contributor

halbmy commented Apr 12, 2024

Dear colleague,
when raising an issue, there is a template that should be filled describing some background, what you did, what you expected and what actually happened. You are copying your code (not even into the part where it should go) along with chinese symbols, no images, but an error message that is clear: You do not have a region with marker 1.

@halbmy
Copy link
Contributor

halbmy commented Apr 12, 2024

You should use only one call of mt.createWorld(). There you can specify layers by the layers argument, see
https://www.pygimli.org/pygimliapi/_generated/pygimli.meshtools.html#pygimli.meshtools.createWorld
and the below linked examples (mostly 2D) using this function.

@zhichaoluan
Copy link
Author

"Sorry, my explanation was not very clear."
image

"I want to generate a 3D layered model and set an anomaly in it to depict the vertical membrane structures in a landfill. During inversion, I want to keep the resistivity of the background model constant and only invert the central vertical membrane area. However, I'm not sure how to set up the 3D layered model, as it seems the 'layer' parameter can't be used in three dimensions. I'm encountering the following error messages in my code."
image

@halbmy
Copy link
Contributor

halbmy commented Apr 13, 2024

Have you tried mt.createWorld(start=[-10, -10, -40], end=[140, 10, 0], layers=[-10, -30])?

Note that you will need a much larger mesh for accurate simulations.

@zhichaoluan
Copy link
Author

"It displays an error like this,"
HV8$UH1{KFA{CS1@_SWWK3F

@halbmy
Copy link
Contributor

halbmy commented Apr 14, 2024

Yes, I am sorry, obviously this feature does not exist yet. I would suggest the following:

  • creating a world as written but without layers (marker 1)
  • adding a cube as a middle layer (marker 2)
  • setting a region marker (marker 3) into the first layer

@zhichaoluan
Copy link
Author

from pygimli.physics import ert
import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
import pygimli.meshtools as mt
world=mt.createWorld(start=[-10, -10,-40], end=[140, 10,0],worldMarker=True,marker=1)
cube3=mt.createCube(size=[128, 0.4, 40], pos=[128/2, 0, -40/2],marker=4)
start_model = world +cube3
pg.show(start_model,alpha=0.3,markers=True)

filename = "cc.shm"
shm = pg.DataContainerERT(filename)
data=ert.load(filename)
data['k']=ert.geometricFactors(data)
data['rhoa']=data['k']*data['u']/data['i']*1000
data.remove(data['rhoa']<0)
data['err'] = ert.estimateError(data,relativeError=0.02)
for s in shm.sensors():
start_model.createNode(s)
for s in shm.sensorPositions():
start_model.createNode(s - [0,0,1e-2/2])
inversion_mesh= mt.createMesh(start_model,quality=1.4)
pg.show(inversion_mesh, markers=True, showMesh=True)
mgr=ert.ERTManager()
mgr.setData(data)
mgr.setMesh(inversion_mesh)
mgr.inv.setRegularization(1, background=True)
mgr.inv.setRegularization(1,fix=200)
mgr.inv.setRegularization(2,startModel=1e4,limits=[100,1e4+20])
mgr.invert(lam=10,zWeight=0.1,verbose=True)
mgr.saveResult()
mgr.showMisfit()

I have created a model that includes an anomaly. In the inversion process, I want to constrain the resistivity of the background to remain constant, optimizing only the resistivity of the anomalous cube. Now, I want to change the uniform background model to a layered model, assigning different resistivities to each layer, and similarly constrain the resistivities of the layered background to remain unchanged during inversion

When I haven't established a layered structure, my forward modeling mesh is correct, as shown below
1713496542290

1713496513170

1713496481079

However, when I try to use a layered model, there is a problem with establishing the mesh; the mesh is not divided correctly
1713496372534

1713496279327
1713496309510

@zhichaoluan
Copy link
Author

cc.zip

@halbmy
Copy link
Contributor

halbmy commented Apr 19, 2024

Please paste code and messages as text instead of screenshots.

In your case the layers (cubes) intersect with the vertical sheet. Whereas such things are automatically resolved in 2D, in 3D this needs to be solved by splitting every cube into two, one on either side of the sheet.

@halbmy
Copy link
Contributor

halbmy commented May 15, 2024

Any news or comments on that?

@kangbow
Copy link

kangbow commented May 28, 2024

from pygimli.physics import ert
import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
import pygimli.meshtools as mt
cube_region_1=mt.createWorld(start=[-50, -10, -31], end=[150, 10, 0],worldMarker=True,marker=1)
cube_region_2=mt.createCube(size=[128, 0.4, 30], pos=[128/2, 0, -30/2], marker=2)
cube_region_3=mt.createCube(size=[200, 9.7, 10], pos=[50, 5.15, -5], marker=3)
cube_region_4=mt.createCube(size=[200, 9.7, 10], pos=[50, -5.15, -5], marker=4)
cube_region_5=mt.createCube(size=[200, 9.7, 10], pos=[50, 5.15, -15], marker=5)
cube_region_6=mt.createCube(size=[200, 9.7, 10], pos=[50, -5.15, -15], marker=6)
cube_region_7=mt.createCube(size=[200, 9.7, 10], pos=[50, 5.15, -25], marker=7)
cube_region_8=mt.createCube(size=[200, 9.7, 10], pos=[50, -5.15, -25], marker=8)
start_model=cube_region_1+cube_region_3+cube_region_2+cube_region_4+cube_region_5+cube_region_6+cube_region_7+cube_region_8
filename = "1D.dat"
shm = pg.DataContainerERT(filename)
data=ert.load(filename)
data['err'] = ert.estimateError(data,relativeError=0.02)
for s in shm.sensors():
start_model.createNode(s)
for s in shm.sensorPositions():
start_model.createNode(s - [0,0,1e-2/2])
inversion_mesh= mt.createMesh(start_model,quality=1.4)
mgr=ert.ERTManager()
mgr.setData(data)
mgr.setMesh(inversion_mesh)
mgr.inv.setRegularization(1, background=True)
mgr.inv.setRegularization(1,fix=100)
mgr.inv.setRegularization(3,fix=150)
mgr.inv.setRegularization(4,fix=150)
mgr.inv.setRegularization(6,fix=200)
mgr.inv.setRegularization(7,fix=250)
mgr.inv.setRegularization(8,fix=250)
mgr.inv.setRegularization(5,fix=200)
mgr.inv.setRegularization(2,startModel=1e4,limits=[100,1e4+20])
mgr.invert(lam=10,zWeight=0.1,verbose=True)
mgr.saveResult()
mgr.showMisfit()

I used the above code to successfully create a grid with three layers, but I ran out of memory on the inversion
1716900042707

1716898859445
But my memory has 180G or so, there should not be insufficient memory problem
1716900214469
How can I solve this problem

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants