Open rowanc1 opened 7 years ago
@rowanc1 this is great! Those two examples are exactly what we need.
That is the model I've been using (meshes behave like lists of primitives) but I've been reconsidering it because it doesn't allow slicing of the meshes, which is very useful for plotting.
My plan was to convert the mesh into something like that indexes like a numpy array:
mesh = TensorMesh([50, 50, 10])
# Get horizontal slice (would be another TensorMesh
layer = mesh[:, :, 0]
# Looping over them would be like looping over arrays
# To get the previous fatiando behavior:
for prism in mesh.ravel():
prism.bounds
This isn't too difficult to do. It would mostly be using slice
objects in __getitem__
and figuring out the geometry of the new mesh.
Would that be useful for SimPEG as well?
There are a couple major differences between the packages, much of the discretize codebase works on vectorization/matrix multiplication to create differential operators. I believe that fatiando is mostly about for loops, and creating (on the fly) representations of the geometric primitives so that calculations can be done on each one?
Our code is designed like that because the meshes and primitives only represent the Earth model. We don't model the fields and data with these meshes as well, like you do in SimPEG because of the numerical PDE approach. We mainly deal with analytical equations for simple mass elements: prisms, polygons, polygonal prisms, spheres, etc.
We also use the primitives to create synthetic data or sometimes on inversions that don't deal with meshes, like estimating magnetization directions.
So our forward modeling is something like:
def gz(x, y, z, model):
result = 0
for element in model:
# calculate analytical solution for element
result += field
return result
So we have exact solutions to approximate geometries. From what I understand, SimPEG does approximate solutions to exact physical property distributions (discretized, of course). This is a nice quirk of fate because it means we have little overlap and can gain a lot from using both libraries together.
There may also be some opportunities to work together on how we deal with physical properties, how these are related to the geometric objects, and how you can abstractly define them (e.g. through mapping classes). This was a big push in SimPEG recently, but that got us to v0.9, still lots of thought to be done on how best to do this.
In Fatiando, they are stored in a props
dictionary within each geometric element. We've come across some examples where this doesn't work because the physical property is somehow tied to the geometry (like in self-demagnetization or depth-varying density). It's also a chore to have to type prism.props['density']
all the time. But I came up with this before I knew of @property
.
I'm curious how you're dealing with these right now in SimPEG. I'm hoping that you've solved this problem already :)
Thanks so much for the response. Lots to think about, I have broken it into a few sections.
Currently for plotting we have opted for a plot function on the Mesh
object. With some modifiers plot_grid
, plot_image
, plot_slice
etc. I am not sure how you are approaching plotting in fatiando.
I see a few things with the __getitem__
sugar. I would like to think forward a bit to unstructured meshes (think classic finite element): array slicing beyond one axis gets weird. I am a bit unclear as to what would be returned by the slicing? Is this a view onto the mesh? Could you do something like:
mesh = TensorMesh([50, 50, 10])
prop = np.array(mesh.ncells) # let's imagine a pep8 future
mesh[:, :, 4].plot(prop) # would show the an x,y image
Your code would likely also work as:
for prism in mesh[10:]: # or just `for prism in mesh`
pass
From my perspective, with the mesh not really being just a collection, the ravel
function on a bare mesh object feels a bit weird (especially for non-logically rectangular meshes, which really only have one axis anyways). I also feel like ravel
would be the default behavior of for x in mesh
.
What is kind of cool about this if that we could do something like only show the cells close to topography or a borehole or something:
cell_nums = compute(mesh, things, topography, borehole)
mesh.plot_grid(faces='k-')
mesh[cell_nums].plot_grid(center='r.', faces='b-')
mesh[cell_nums].plot(prop, opacity=0.8) # notice that we did not do `prop[cell_nums]` here
Interesting, composable and pretty clean. The plotting inside discretize
needs some love regardless and moving some things around to make it a bit more clear/approachable.
Regarding mag/grav and analytical functions from the meshes/geometries, you should touch base with @fourndo, who has done a bunch of work in the SimPEG code base for both analytic and differential versions of both of these physics. Each has their benefits - better to have both, and just let them work exactly the same way. :)
I believe that he is currently doing some Magnetization Vector Inversions, as well as joint inversions etc. @fourndo moves too fast for me to keep up.
We addressed parts of it in simpeg/simpeg#444, and have gone through at least three ways of doing this as our thinking has evolved over the last four years! There is a video that actually walks through this issue, it is a bit dense.
https://www.youtube.com/watch?v=VtFC7jQVWaM&feature=youtu.be&t=5m55s
We are looking towards multi-physics joint inversions, for example, think a flow problem monitored by an EM problem, but also having some in situ saturation measurements, with a distributed temperature sensor or something (this is now!). And you want to use all the data and invert for hydraulic conductivity, or saturation, or electrical conductivity, or both while regularizing with a cross-gradient term, or invert for the n
parameter in Archie's equation. Or you want to run the whole thing in the forward mode, and not care about any inversion crap. Or you want to choose a model
to be some sort of low dimensional parameterization, like a boundary between two layers or a Gaussian contaminant plume.
So a few constraints:
Because we think a lot in matrix equations, these properties end up being arrays with their index relating to the cell numbering (if attached to a mesh). They are also (currently!) "dumb" on their own. They don't really know if they live on a mesh, can't plot themselves, and are just numpy arrays.
We have also named these properties on the problem/physics object instead of on the mesh. For example, the e-formulation of Maxwell's equations Maxwell_e
has a concept of sigma
(electrical conductivity) - the mesh has no clue about those things.
mesh
- holds operators (differential, averaging), can yield geometric primitives (only the TreeMesh
at the moment actually does this)geometric_primitive
- don't really have this implemented everywhere, helpful in plotting & integral physicsprop
- (physical) property that is at some location on the mesh/geometric-primitive (e.g. center, nodes, faces), may not be attached to the concept of a mesh, makes more sense on a Problem/Physics
plotting
- any combination of the above three.Bit of a brain dump, hope it helps rather than scares you away. :)
We are just getting into experiments with v4 of this concept, I am relatively optimistic at this stage.
Looking forward to hearing more.
@rowanc1 I haven't forgotten about this, just taking some time to think about it more deeply :smile:
@rowanc1 I've been giving this a lot of thought lately. I'm not happy with the way Fatiando is handling all this but I'm not sure how I could fit my current mentality into the separations of mesh/problem/physics you have in SimPEG.
I've been trying to make all my assumptions explicit so that I can start to understand how this all fits together. For that, I'm (very, very slowly :) stripping out fatiando.mesher
into its own package fatiando/geometric and rethinking how I want to handle physical properties.
I'm still experimenting with implementations and ideas but I'll write about how this thing works when I arrive at something that solves my current problems. I think this will be a good starting point for a later merge with discretize.
Regarding plotting, all I can say is: default to matplotlib for everything. Even if it's not pretty, everyone can get it installed. Unless you really want to get a lot of complaints about Mayavi installation 😉
I am leaving this open for a general open discussion around the idea of bringing some aspects of the fatiando meshing/gridding software to either make use of or replace aspects of this package. With the goal of de-duplicating small parts of the geo-scientific stack and bringing in some common dependencies.
The
SquareMesh
for example is quite similar to theTensorMesh
in 2D: http://www.fatiando.org/api/mesher.html#fatiando.mesher.SquareMeshThe 3D tensor mesh and octree meshes could generate
PrismMesh
es: http://www.fatiando.org/api/mesher.html#fatiando.mesher.PrismMeshIt looks like there would need to be some work in looping over cells to generate geometric primitives. This is actually exactly what we are doing in the treemesh code (this code is dense). https://github.com/simpeg/discretize/blob/master/discretize/TreeMesh.py#L2282
But basically what I could see is:
There are a couple major differences between the packages, much of the discretize codebase works on vectorization/matrix multiplication to create differential operators. I believe that fatiando is mostly about for loops, and creating (on the fly) representations of the geometric primitives so that calculations can be done on each one?
There may also be some opportunities to work together on how we deal with physical properties, how these are related to the geometric objects, and how you can abstractly define them (e.g. through mapping classes). This was a big push in SimPEG recently, but that got us to v0.9, still lots of thought to be done on how best to do this.
Again, I see this as a way to bring our two communities a bit closer together, as there are significant overlaps: grav/mag, srtomo, inv framework, utilities like data fetching, etc. But one step at a time. :)
@leouieda has noted some of these things in fatiando/fatiando#278 and also pointed out a few opportunities to improve the
discretize
package (e.g. pep8) #39, which we should seriously consider.