Block model embedded in 1D VTI background

[prisae]

1D Model

You can find the 1D model in the notebook notebooks/BlockModel.ipynb.

Please have a look and give your feedback (@lheagy, @Moouuzi, @ocastilloreyes). If everybody is fine with it then we can fix this, so everyone can try to model it with their codes.

I used a rather coarse grid now for testing, so precision is not overwhelming right now (but with roughly 2 minutes reasonably fast to calculate). I might refine this for emg3d once we agree on the model.

Questions:

• How do you need the data? Is this notebook sufficient?
• Should we make it VTI? Can all the codes handle VTI resistivity?
• Currently x-directed dipole, 200 m long. Could make that an arbitrary rotated dipole - can all codes handle arbitrary rotated dipoles?

Model and survey lines

1D background, compared to empymod

3D block model

(Related to #4.)

Now I will try to get the Marlim R3D model working...

[Moouuzi]

I just built the mesh as it is now: In general, very thin layers over 20 x 20 km^2 horizontal extent with good-quality tetrahedra demand a lot of elements. I'm not saying that I cannot calculate this, but as this should be somewhat easily reproducible from other users, also with unstructured codes, I'm either for increasing the thickness of the 5 Ohmm layer to 100m or maybe removing it at all. If it is still reasonable, I would prefer even even greater thicknesses, like 200 and 100-150 or 300 and 200 m for the first 2 layers below the seawater. The following numbers are still without receiver refinement and anomalies, but this requires only a "small" number or additional nodes...

Now: 500 k nodes --> p1: 10 min and ~ 500 GB RAM or p2: 2h and 3 TB RAM (absolute maximum) 100 m thick 5 Ohmm layer: 150 k nodes --> p1: 1-2 min, p2 10-20 min no 5 Ohmm layer: 30 k nodes -_> p1 <10 s, p2 1 min

Aside from that, arbitrary sources/anisotropy is no issue for me.

Another comment: The misfit of SINGLE components (x, y, z) in oblique directions or on oblique observation lines is about 2-3 times higher than for orthogonal setups. Of Course this makes sense if some components are zero, but it also depends on the strucuture of mesh design (direction of dof and interpolation around the receivers). Therefore, I would prefer an oblique transmitter, if possible for everyone.

[prisae]

@Moouuzi thanks for the feedback. Making this top layer much thicker sounds good to me. I will wait for the feedback of @lheagy and @ocastilloreyes , before making changes though.

So plans for now, if the others agree:

• Increase thicknesses of the top layers
• Implement slight VTI anisotropy (I know SimPEG can handle it, what about PETGEM?)
• Rotated source, what about 20° azimuth and 10° dip? (Again, I know SimPEG can handle it, what about PETGEM?)

Working on the MarlimR3D model now, will soon post some images on #10 , it looks quite interesting!

[ocastilloreyes]

My comments are:

• I agree to increase the thickness of the top layers.

• Yes, PETGEM supports VTI. For this case, I will need two conductivity values per element (sigmaX == sigmaY, sigmaZ)

• Rotated sources are not supported. The current version only considers x-directed, y-directed or z-directed dipoles (I had planned to include this point in the list of features to be developed in the future, my wish list..). However, I think it is quite simple to implement the source rotation. I will review it and define my reference/rotation system. At the end we will discuss the compatibility of rotation systems. What do you think?

[Moouuzi]

@ Octavio

If you have all x-y-z directed Sources implemented, it's quite easy to rotate them. It's like projecting on the reference element for intgration. i just asked my colleague Nico, he is doing it as you already thought about, you just need a simple rotation matrix:

1. just horizontal rotation (x-y): y> rotate all corrdinates -> calculate for x-directed HED -> rotate fields back to match the original orientation

2. full rotation (z): calculate horizontal component as in 1.-> calculate vertical component with z-directed HED -> sum up resulting source vector.

[ocastilloreyes]

@Moouuzi great idea! Thank you. I will check/implement this approach. If it works, yes, PETGEM will "support" rotated sources.

On the other hand, I have some questions that can be applied to the 1D and 3D model.

• How will we manage the mesh creation?
• For FEM codes, What mesh format will we use? Is it possible to use the same mesh for different FEM codes?

I think we should try to take advantage of the common points of our codes. In this way we can focus on the particularities of each one. This will save us time / work. But of course, I understand that this could be complicated to achieve. What do you think?

[Moouuzi]

• using a common mesh format should be easy, I'm currently working with TetGen, which exports .mesh (Medit) or vtk, and convert them to fenics-conformal XML or h5 meshes.

• I already offered Dieter to take care of the 3D tetrahedral mesh generation, of course you could also use my meshes and we don't have to spend the effort twice, especially since you have other duties now:-)

• big pro: saves time, differences in codes are not related to geometry-differences in our meshes

• big contra: If the paper gives an overview about open-source codes, but users are completely left alone with mesh generation for a complex model, they are left alone with the most challenging, time-consuming and accuracy-controlling part. I think we all agree that pressing the "run and interpolation buttons afterwards" is mostly automated in all of our codes.

My suggestion to start: If I manage to build an approximated mesh from the 3D Marlim example with my mesh-tools as they are, I would use this example to write a tutorial for creating real-world related meshes anyway and we can add a jupyter Notebook to the publication, which enables users to use this mesh with either PETGEM or custEM to recompute our results. If you have alternative ideas for creating the mesh in the meantime, even better.

[prisae]

Agree regarding mesh. A part of it could be improved with GemPy, we tried that already on an example, and it works nice.

And yes, rotated sources are fairly easy, so you should not encounter problems. At least as long as the entire source is in one medium and not crossing boundaries. Otherwise it only becomes slightly more difficult, as we have to take care of the different medium parameters the source touches.

Cool, slowly things get moving. Let's see how it goes. Otherwise it might be worth having a call once to discuss a few things.

[prisae]

OK, I updated the 1D model with your inputs. Three changes:

• Now one thick top-layer of 250 m.
• Below that, the background has a slight VTI anisotropy.
• Source is a slightly rotated, finite length dipole.

[lheagy]

This is looking solid! Sorry I have been a bit MIA on these conversations. I think including VTI is interesting. We certainly can do rotated sources, but I am wondering if it might be simpler in the discussion and description of responses if the dipole is axes-aligned. I think we may want to point out a few features in the data, and rotating the source does slightly complicate things.

If this is something you feel is important to highlight though, that is fine by me!

[lheagy]

Also, just a quick clarification - should we re-name this issue to "block model"? The model is 3D, and I presume you are running the simulation with 3D codes, right?

[prisae]

When I tried I also thought it might be better to have just an inline source, as in the end a rotated source is only a linear combination of aligned sources anyway. But @Moouuzi mentioned that it might be better to have a rotated source:

Another comment: The misfit of SINGLE components (x, y, z) in oblique directions or on oblique observation lines is about 2-3 times higher than for orthogonal setups. Of Course this makes sense if some components are zero, but it also depends on the strucuture of mesh design (direction of dof and interpolation around the receivers). Therefore, I would prefer an oblique transmitter, if possible for everyone.

Maybe I understood him wrong though. @Moouuzi , would you prefer a rotated dipole or not?

[Moouuzi]

A rotated source is not realy necessary, it just makes the setup "realy 3D" in terms of no non-zero components anywhere.

On the other hand, Lindsey is right, if you want to point certain features, orthogonal alignement of anomalies, tx and rx is of advantage. In addtion, the model should not be too complicated to be remodeled. Meanwhile, I think an orthogonal setup is very adequate - this is the first example and the Marlim model will be completely 3D anyway.

[prisae]

Okay. I updated the model to have an x-directed, 200 m long dipole.

=> Any other concerns, or can we call this pre-final? (I guess only once everybody calculated it and we compared it and everyone is happy with it we can call it final.)

[ocastilloreyes]

I am online again! After my paternity leave I am somewhat lost. I hope to be in tune again ...

Today I implemented the source rotation, I validated it and everything is ok. My idea is to work on this modeling next week (the rest of this week I will dedicate to finish some pending tasks). Right now I have a question about the mesh for this block model: @Moouuzi Do you have any mesh for this case or should I build a mesh?

[Moouuzi]

Hi Octavio,

I hope you had a nice time and everyone is healthy! I still need to resolve a small issue for adding the "intersecting" blocks to my poly-file for TetGen, but that is fixed very soon. I'm fine if you want to use my mesh or the poly file. The only question would be that readers might be interested to see how such a mesh can be generaten with PETGEM utilities, as this is not the complicated 3D world example. Anyway, I will send you the poly-file by end of the week, and also the one for the marlim mesh.

Cheers, Raphael

PS: Thank you for forwarding the travel grant call, I just had no time to deal with it until now.

[Moouuzi]

Finally, I also managed to build the blocky mesh. Intersecting blocks in combination with a layer intersection are quite challenging for building a confrom TetGen input. It's not a perfect solution, as I need to devide the block into sub-blocks that have identical sub-surfaces for the autmated mesh generation, but that does not affect the quality.

I will run the computations in the next days. Do we already have a directory for uploading data for the comparison? Or should we send our results at the Rx locations to you Dieter and you are calculating misfits and generating the plots?

[prisae]

Not quite sure what you mean with intersecting blocks - none of them is overlapping the other, the are just touching.

Well, here, there are not too many receivers, right? So I guess you could just upload it here on GitHub. Just create I directory data, and I will put mine too and I'll compare a first comparison plot next week with it!

Looks awesome, looking forward to the comparison! Just seeing the models coming together I think already just that will be very useful to the community!

[Moouuzi]

Sorry for the misunderstanding, of course, the blocks are just touching each other. For me, "touching" and "intersecting" interfaces require an identical handling, as both cases results in "intersection" errors from TetGen, independent of the fact if faces are touching or intersecting each other.

[ocastilloreyes]

A few minutes ago I created a new branch for petgem results. Please look at the results with p1. The interpolated electric field results are located in the data/petgem_results_blockModel.

The file format is the same as @Moouuzi has used (nx3 arrays). I'm interested:

• Verify that petgem modeling parameters are correct

• The negative effect of boundaries. I think that the computational domain should extend in the horizontal direction, but this can be discussed.

The computation times are ~120 seconds with 24 cores.

I think we are about to close this model. To cross your fingers!

[prisae]

Great, I reserved this afternoon for that, so I will compare then and hopefully they are already close enough, otherwise we iterate once over it! I'll let you all know!

[prisae]

Awesome, thanks @Moouuzi and @ocastilloreyes , we have our first comparison plots, and they look already awesome! (Only E-field at the moment, I should add H-field at some point.) Now just SimPEG is missing (@lheagy).

I merged both your results, and I created a new branch, https://github.com/prisae/3d-csem-open-source-landscape/tree/blockmodel-results => I hope you all have write access to that, so you could make changes straight in there for the Block Model.

We'll also have to make a table. The table has to contain:

• how many cells
• runtime
• how many threads
• On what type of machine run

=> However, I don't think it makes sense to do that now. Mostl likely we'll run them again before submitting, once we have a draft of everything. So we'll collect the final results when we are there.

Once we have a working example for Marlim too and we know that it is all going to work out nicely I'll shall draft a paper, and we iterate over the paper before we spent more time on (time-consuming) modelling...

[prisae]

For the comparison, see: https://github.com/prisae/3d-csem-open-source-landscape/blob/blockmodel-results/notebooks/BlockModel-Comparison.ipynb

[Moouuzi]

Nice first results!

I think, I uploaded the p2 results. The horizontal extent for the p2 mesh is very small, also not too large for the p1 mesh. I wanted to keep the computation times low for the first run, but I agree with Octavio that the extents ars not sufficient for both meshes, that's why the boundary artifacts. So either we switch to p1 or even lower Quality for p2. I will make some tests to see which extent is somewhat sufficient.

We could also think about analytic boundary conditions in the 1D case, but I don't have them available for anisotropy. I would require a wrapper-function for empymod, that takes the coordinates of the boundary dof and calculates the semi-analytic field solutions on these points. With that, the mesh just needs to be as large as the receiver extent. Anyway, this does not help at all with the more complicated 3D models, so I think it's not worth the effort.

[ocastilloreyes]

Yeah! the results are nice...

@Moouuzi I usually need to extend the domain at least four times the smallest skin-depth value. In marine modeling, this minimum value commonly occurs in water (for this case I have obtained min_skin_depth = 273.86m). However, I have only verified this strategy for isotropic models and I have no reference on how it would work in cases with anisotropy. In addition, this depends on the numerical scheme. Therefore, I think this is a good opportunity to try it and understand it!

I also think it is not worth the effort to use/implement analytic boundary conditions... In any case, I'm glad to know that this is moving on. Thank you all!

[prisae]

It is probably not worth it, no. Particularly as we want to focus on the 3D modellers. Depending on space we might not even include the pure 1D results in the publication, but refer to it and having it in the notebooks/repo.

As info: If you want to have really nice results then you would probably need 6 skin depths. And, further more, take the skind-depth into account for air, not just with upwards pointing z, but also horizontal (x and y), because the airwave gets bounced back and hit on the sides as well. I just finished some testing on this. Usually the airwave is only taken into account vertically, which is wrong.

For the relatively precise emg3d results these are the rough boundaries:

• x: -106 km - 106 km
• y: -102 km - 102 km
• z: -34 km - 144 km (positive z upwards)

This is certainly extreme. But it helped to keep relative error below 2% (except by zero crossings).

[Moouuzi]

@Dieter

thanks a lot for these additional information! Usually, I also use ~100 km extent for most modelings, but in this layered-earth case, this is way too big for a good-quality tetrahedral discretization! I will also try to find out if a halfspace-boundary mesh helps in this case

[prisae]

There is something else that might help and speed up: You just need the air-subsurface interface for +/- 100 km, there wont come anything back from the subsurface. So sufficiently far out you can simply get rid of the layered subsurface model, and only have one interface (air above and the resistivity of the top layer below). That would be interesting to see if it works, but I think so.

Either way (this approximation, analytical approximation, etc), we'll have to make sure to document them properly.

[Moouuzi]

That's exactly what I'm doing for a long time and I wanted to express this approach with "find out if a halfspace-boundary-mesh helps" :-) It works in land-based CSEM modelings quite well.

[prisae]

Ah Ok. Yes, I would expect it works. In this case just air above water. If you plot entire 2D slices of the field, you see that there is only energy in the air (far out). So it goes in the air to the boundary, then back, and then basically down to the receiver. So nothing will come back horizontally through the water, and even less through the subsurface. Just src->up->through-air-to-boundary->through-air-back-above-receiver->down-rec.

[Moouuzi]

Dear all,

I optimized the mesh design for the 1d case - in my branch, there are new p1 and p2 results witch require also about 2-3 minutes computation times, but have way less artifacts due to appending the halfspace-like boundary mesh. I will also update the meshes for the blocky model and forward the new mesh-files to Octavio. One thing I can already say: it will be hard to reach the 1% error level with first order polynomials.

I wish you all a relaxed vacation time! Raphael

[prisae]

Cool, great!

due to appending the halfspace-like boundary mesh

Help me to remember this, as this is useful/important info we'll have to describe in the paper.

hard to reach the 1% error level

I guess we don't have to be too strict on that, and 2-3 % is probably fine. Again, not worth bothering too much in detail yet. Once everything falls into place we'll see if it is good enough or not.

[Moouuzi]

Hey, I finally updated also the block-model results. @ Octavio: mesh generation files will be sent in a few minutes @ Dieter: could you copy the results from the custem_input branch to the master branch and update the plots in the notebook? I would call this my final version.

[prisae]

Thanks! I'll do it tomorrow.

[prisae]

Hm, I see that your branch 'custem_input' is a mix of block-model and 3D results etc, which makes it difficult to merge parts of it...

I will merge now the entire branch. Please delete the branch, and then down the road it will be easier if you make a branch for every different thing, so we can easily create PRs and merge particular things.

[prisae]

So I merged all open branches into master. There are two branches which I did not delete, I let you (@Moouuzi , @ocastilloreyes ) do that, to not delete something by accident:

• custem_input
• petgem_input

So everything is in master. Please update your repo, and then create new branches for each feature/result/text. Thanks!

[prisae]

(I'll write the rest on Slack.)

[ocastilloreyes]

Hi, I have finally updated the results for the block-model (my final version)

@Dieter could you do the merge please?

[prisae]

I merged it into master. However, when I run https://github.com/prisae/3d-csem-open-source-landscape/blob/master/notebooks/BlockModel-Comparison.ipynb your data looks very different, so I probably do something wrong. Can you have a look at the notebook and tell me how to adjust the PETGEM block?

# BACKGROUND
ptg0 = np.load('../data/petgem_results_blockModel/y_0_E_bm.npy')
# ptgp = np.load('../data/petgem_results_blockModel/y_3000_E_bm.npy')
ptgm = np.load('../data/petgem_results_blockModel/y_-3000_E_bm.npy')
ptg_bg = np.vstack([ptgm[:, 0], ptg0[:, 0]]).T
# Take conjugate, PETGEM uses -iwt, as oposed to emg3d/custEM which have iwt.
ptg_bg = ptg_bg.conj()

# 3D
ptg0 = np.load('../data/petgem_results_blockModel/y_0_E_le.npy')
ptgp = np.zeros((101, 3), dtype=complex)                                          # OCTAVIO, y_3000_E_le has one entry less?
ptgp[:-1, :] = np.load('../data/petgem_results_blockModel/y_3000_E_le.npy')
ptgm = np.load('../data/petgem_results_blockModel/y_-3000_E_le.npy')
ptg = np.vstack([ptgm[:, 0], ptg0[:, 0], ptgp[:, 0]]).T
# Take conjugate, PETGEM uses -iwt, as oposed to emg3d/custEM which have iwt.
ptg = ptg.conj()

One adjustement is obviously to append _p1 or _p2, but other than that?

[ocastilloreyes]

@Dieter Sorry, it was my mistake. The data/names were reversed. I have already solved it and verified it in your notebook. The updated/correct data is in my branch.

Thanks for finding the error!

[prisae]

Great, I just updated https://github.com/prisae/3d-csem-open-source-landscape/blob/master/notebooks/BlockModel-Comparison.ipynb on master, looks great now! For the results with blocks PETGEM and custEM are basically identical, emg3d is a bit off still. But these are still my first quick try results when I created the model, so I will have a more in-depth go at it to match you results.

@ocastilloreyes: Did you see my comment in the notebook: "# OCTAVIO, y_3000_E_le has one entry less?"

@ocastilloreyes @Moouuzi : Can you post here (or in a file on your branches) the computational information:

• Total CPU time
• Total max RAM
• How many cores/threads
• What hardware
• Other info that is relevant for FE methods, which I don't know enough (nuber of cells or nodes or edges or ...?)

@lheagy How is it going with SimPEG?

[ocastilloreyes]

@prisae sorry .. I have already fixed the error with the y_3000_E_le entry (correct data is in my branch).

About computational information, of course, during the weekend I will upload a file to my branch

[prisae]

No problem. Included the new data, all good now.