gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
374 stars 137 forks source link

Robin BC in [Geoelectrics in 2.5D] #269

Closed HenryWHR closed 3 years ago

HenryWHR commented 3 years ago

Hello everyone

I was trying to understand the mixed boundary condition in resistivity modeling with some examples and I found the Geoelectrical in 2.5D is a good one. However, I have some questions.

  1. In the example, to determine the alpha for the robinBC parameter, sigma is also involved. With some tests, I can see this sigma here matters. alpha = sigma * k * ((r1.dot(n)) / r1A * pg.math.besselK1(r1A * k) + (r2.dot(n)) / r2A * pg.math.besselK1(r2A * k)) / \ (pg.math.besselK0(r1A * k) + pg.math.besselK0(r2A * k)) However, the formula for this parameter presented in Dey & Morrison, 1979 and some other references seems model-independent, which confused me.

  2. I wonder how the robin BC is applied for heterogeneous models. In my understanding, the potential values at the surrounding boundaries are derived from the homogeneous analytic solution a prior, even for the heterogeneous model. I am not sure if this is true. But, if so, which homo model should be used to obtain the boundary value.

I appreciate your help with my questions. Many thanks!

Best regards

Henry

halbmy commented 3 years ago

Hello, generally, I think that mixed boundary conditions make only sense for homogeneous models, where the behaviour (i.e. the angle the current flow toughes the boundary) is known. For heterogeneous models like layers, one would rather choose inhomogeneous Dirichlet conditions.

The matter is easy in 3D, where the mixed BC basically represent the mentioned angle. Therefore, for 3D the alpha value is basically a cosine (dot product) between source vector and normal. (Note that we use a median source position to avoid reassembling any matrix, the boundaries are far away from the region of interest and we model electrodes separately anyway.)

However, in 2D we have the wavenumber term with the additional u so that sigma determined the decomposition and therefore the boundary conditions. I've looked 2-3 references but it is always the same for 2D, if you find any without sigma in 2D, please tell me.

HenryWHR commented 3 years ago

Hello, Thank you very much for answering my question. However, I still don't really understand. There are some references below where I didn't see the alpha has a sigma. Those equations are almost the same (some for surface electrodes only). I guess I missed some math here?

Page 113 in Ref [1], Page 47-48, equation (3.8), (3.13), (3.14) in Ref [2], Equation (10) in Ref [3], Equation (2.5c) in Ref [4].

[1]. Dey, A. and Morrison, H.F., 1979. Resistivity modelling for arbitrarily shaped two‐dimensional structures. Geophysical Prospecting, 27(1), pp.106-136. [2]. Kemna, A. 2000. Tomographic inversion of complex resistivity. Ruhr-Universität Bochum. [3]. Pan, K. and Tang, J., 2014. 2.5-D and 3-D DC resistivity modelling using an extrapolation cascadic multigrid method. Geophysical Journal International, 197(3), pp.1459-1470. [4]. Yuan, Y., Qiang, J., Tang, J., Ren, Z. and Xiao, X., 2016. 2.5 D direct‐current resistivity forward modelling and inversion by finite‐element–infinite‐element coupled method. Geophysical Prospecting, 64(3), pp.767-779.

Thanks again for helping me through this.

Best regards

Henry

halbmy commented 3 years ago

I totally agree and I was wrong before. The wavenumber decomposition cannot change the fact that the ratio between the potential and its derivative is purely geometrical and thus not depending on sigma. In Kemna (2000) all is very clear.

However, the implementation in dcfemmodelling.cpp is different. Maybe @carsten-forty2 can shed some light on it.

halbmy commented 3 years ago

Well, the boundary condition, i.e. the ratio between u and du/dn is not depending on sigma. But the entries in the stiffness matrix are conductances, so the change in the A matrix needs to contain sigma as well (if using a += notation, it would be different using a *= notation).

halbmy commented 3 years ago

All is well, if you look at (3.22) in Kemna (2000): There is a sigma for the elements and also a sigma (in beta) for the boundaries so it does not change the relation between u and du/dn, but of course sigma goes into the stiffness matrix.

HenryWHR commented 3 years ago

Thank you very much!

now, I think I understand. So, if I understand correctly, there should also be a sigma in 3D case as this 2.5D.

halbmy commented 3 years ago

If you just talk about the ratio between the potential and its normal derivative, there is no sigma, neither in 2D nor 3D.

However, when implementing the BC, one needs to change the matrix entries of A, and they are conductances and thus include sigma, both in 2D and 3D. It is a matter of implementation if it is in the according function or not.

HenryWHR commented 3 years ago

Yes, thank you very much for your help. I think I understand now.

Best regards

Henry