Fatiando

[rowanc1]

Fatiando

[leouieda]

@rowanc1 I haven't forgotten about this, just taking some time to think about it more deeply 😄

[leouieda]

@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 :)

[rowanc1]

Thanks so much for the response. Lots to think about, I have broken it into a few sections.

---

Plotting & Slicing

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.

Mag/Grav

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.

Properties

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:

• need clear separation between the model and the properties (connected by some sort of mapping and associated derivative)
• properties can either look to the model, or know what they are already
• have defaults when they make sense, yell when you screw up etc.
• be able to calculate derivatives if there is model dependence
• may not be tied to a mesh (could be a parameter in a source term)
• be simple to use

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.

Concepts

• 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 physics
• prop - (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.

---

Overview

• I think there is a space to have the geometric primitive concepts better defined
• I am coming around to the idea of slicing into the mesh. I would like comments from @lheagy
• This is an opportunity to do a better job at plotting, and take some of your learnings/thoughts

Looking forward to hearing more.

[leouieda]

@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 😉