ruboerner/femt2d

FE simulation of 2D Magnetotellurics in MATLAB

FEMT2D

This repository contains the MATLAB toolbox FEMT2D.

It provides a discretization of Helmholtz type PDE's in two dimensions using first or second order Lagrange elements on triangular meshes.

Requirements

Besides a basic MATLAB installation, a mesh generator is required for setting up the triangulation. Currently, only triangular meshes generated by the Triangle library (references see below) are supported.

To make the creation of 2-D triangular meshes more user-friendly, the installation of PyGIMLi is strongly recommended.

FEMT2D has been tested with MATLAB R2020a. Earlier releases work well also. Also note that compatibility with Octave has not been tested.

To activate rendering of $\LaTeX$ expressions within this README.md, download and activate this Chrome extension.

References

• FEMT2D: FE simulation of 2D Magnetotellurics in MATLAB. https://doi.org/10.5281/zenodo.3955407
• Jonathan Richard Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, in "Applied Computational Geometry: Towards Geometric Engineering" (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of Lecture Notes in Computer Science, pages 203-222, Springer-Verlag, Berlin, May 1996.
• Jonathan Richard Shewchuk, Delaunay Refinement Algorithms for Triangular Mesh Generation, Computational Geometry: Theory and Applications 22(1-3):21-74, May 2002.
• Franke, A., Börner, R. U., & Spitzer, K. (2007). Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography. Geophysical Journal International, 171(1), 71-86.
• Börner, R. U. (2010). Numerical modelling in geo-electromagnetics: advances and challenges. Surveys in Geophysics, 31(2), 225-245.
• Rücker, C., Günther, T., Wagner, F.M., 2017. pyGIMLi: An open-source library for modelling and inversion in geophysics, Computers and Geosciences, 109, 106-123, doi: 10.1016/j.cageo.2017.07.011.

Installation

Clone this repository.

git clone https://github.com/ruboerner/FEMT2D.git

Then enter the newly created folder FEMT2D:

cd FEMT2D

Quick start

For a quick demonstration, simply call

>> demo.driverWeaver

which loads a triangular mesh, and evaluates the response of two quarter-spaces for a given period at predefined observation points. A few typical plots, such as profiles of apparent resistivities, phases, and magnetic transfer functions will be generated.

Folder structure

Inside the FEMT2D folder, there are a few subfolders. The folders with leading +, e.g., +fe, contain packages and have their own namespace.

FEMT2D
├── +demo
├── +fe
├── +mesh
├── +mt
├── +plot
├── +tools
├── Notebooks
├── meshes

From the top-level folder FEMT2D you can start a demonstration by typing, e.g.,

demo.driverCOMMEMI_2D_0

which will calculate the model responses of model 2D-0 from the COMMEMI model suite.

The folder Notebooks contains example Jupyter Notebooks which introduce the user into the process of mesh generation using PyGIMLi.

Each of these Notebooks creates a particular set of files in the folder meshes. The meshes can later be read by routines which reside in the +mesh folder.

Details of the toolbox

After reading in a Triangle triangulation file, the mesh topology is stored within a MATLAB struct.

mesh = mesh.getMesh('filename', 'meshes/commemi2d0.1', 'format', 'triangle');

'filename' corresponds to the basename of a mesh file in the meshes folder. Note that the actual files have the extensions .poly, .ele, and .node. If the files in the meshes folder are, e.g.,

meshes
├── commemi2d0.1.ele
├── commemi2d0.1.node
├── commemi2d0.1.poly

then the corresponding filename in the parameter list of mesh.getMesh() would be commemi2d0.1.

The mesh structure

The mesh structure contains the following fields:

field name	dimension	description
node	2 x np	node coordinates
tri2node	3 x nt	table of triangles defined by their indices as listed in the array node
tri2subdomain	1 x nt	array of subdomain number each triangle belongs to
marker	1 x nt	array of attributes associated with nodes
edge2node	2 x ne	table of edges and their node indices
tri2edge	3 x nt	table of edges and their associated triangle index
nt	scalar	number of triangles $n_t$
np	scalar	number of nodes $n_p$
ne	scalar	number of edges $n_e$
bdEdges	vector	indices of boundary edges
bdNodes	vector	indices of boundary nodes
Bk	2 x 2 x nt	array of matrices defining the affine map $B_k$ from reference element to each of the $k=1,\dots,n_t$ domain elements
Binv	2 x 2 x nt	array of matrices $B_k^{-1}$
detBk	nt x 1	array of determinants of $B_k$, $\text{det}(B_k)$
area	nt x 1	array of all $k=1,\dots,n_t$ triangular element surface areas, i.e., area = 0.5 * detBk

|

The fem structure

The finite element method associates the triangulation with specific basis functions defined on the elements and identifies the degrees of freedom (DOFS).

freq = 0.01;
sigma = [1e-1, 1e-1, 1e-2, 1e-2, 1e-14];
mu = ones(size(sigma));

fem = fe.FEMproblem('mesh', mesh, ...
    'elementtype', 'Lagrange', ...
    'order', 2, 'dimension', 2, ...
    'sigma', sigma, 'mu', mu, ...
    'polarization', 'both', ...
    'frequency', freq, ...
    'verbose', true);

As one can conclude from the above example, the mesh consists of five subdomains each of which is associated with a unique parameter (here sigma denotes the electrical conductivity $\sigma$ given in S/m and mu is the relative magnetic permeability $\mu_r$, such that $\mu = \mu_r \mu_0$ with $\mu_0 = 4 \pi \cdot10^{-7}$ Vs/Am).

In the given example, the solutions of the Helmholtz equations (for both magnetotelluric polarizations) are desired at a frequency freq of 0.01 Hz. The polynomial order order of the Lagrange elements is 2, i.e., we use quadratic (second order) finite element basis functions. The numerical accuracy is always much better when second-order elements are employed.

After calling FEMproblem, the fem structure contains the following fields:

field name	dimension	description	default
mesh	struct	copy of mesh struct	-
order	scalar	1 or 2 for linear or quadratic elements, resp.	1
dimension	scalar	spatial dimension (currently only 2D is implemented)	2
app	struct	contains application (PDE) specific parameters, e.g., the frequency frequency for the MT case	-
polarization	string	indicates whether E- ('epol'), H- ('hpol') or both ('both') polarizations are desired	'both'
elementtype	type of finite elements	'Lagrange'
sigma	nt x 1	elementwise electrical conductivities $\sigma$ in S/m	-
mu	nt x 1	elementwise relative permeabilities $\mu_r$	-
stripair	scalar	indicator that controls whether the Air halfspace has to be excluded (e.g., for the MT for H-polarization)	-
istrip	scalar	subdomain number associated with the Air halfspace that has to be stripped from mesh for MT	-
coorddofs	2 x ndofs	coordinates of the DOFS	-
elem2dofs	(order * 3) x nt	array of elements and their associated DOFS	-
ndofs	scalar	number of degrees of freedom (DOFS)	-
bddofs	vector	boundary DOF indices	-
Compl	ndofs x nbddofs	mapping from full to boundary DOFS	-
Null	ndofs x (ndofs - nbddofs)	mapping from full to inner DOFS	-
EnL	scalar	electric field at $z=0$ for the left boundary	-
EnR	scalar	electric field at $z=0$ for the right boundary	-
HnL	scalar	magnetic field at $z=0$ for the left boundary	-
HnR	scalar	magnetic field at $z=0$ for the right boundary	-

The assembly of the FE matrices can be carried out using a call to FEMassemble:

fem = fe.FEMassemble(fem, 'output', 'matrices', 'verbose', true);

After return, the fem struct contains the following additional fields:

field name	dimension	description
Ke	ndofs x ndofs	E-pol stiffness matrix (sparse)
Me	ndofs x ndofs	E-pol mass matrix (sparse)
Kh	ndofs x ndofs	H-pol stiffness matrix (sparse)
Mh	ndofs x ndofs	H-pol mass matrix (sparse)

The observation operators

Usually the solution (and its gradient) is only required at a few points within the computational domain. To this end we setup observation operators for the desired polarizations and field components, such that $\mathbf Q : \mathbb C^N \to \mathbb C^n$ maps from an $N$-dimensional discrete solution (or components of its gradient) to an $n$-dimensional vector of spatial observations. Usually, $N \gg n$ and $\mathbf Q$ is sparse. The observation operators are independent of the frequency.

The observation operators can be constructed using the call

yobs = tools.asRow(-99000:1000:99000);
obs = [yobs; 0.01 + zeros(size(xobs))];
nobs = size(obs, 2);
fem = fe.getQ(fem, obs);

The fem structure now contains an additional field Q which in turn contains the discrete observation operators.

It is always assumed that the strike direction is aligned with the $x$-axis of a right-handed coordinate system. More precisely, the $x$-axis points into the plane, the $y$-axis point to the right, and the $z$-axis points downward.

field name	dimension	description
Qe	nobs x ndofs	evaluates $E_x$ for E-pol
QeHy	nobs x ndofs	evaluates $H_y$ for E-pol
QeHz	nobs x ndofs	evaluates $H_y$ for E-pol

To enable the plot of numerical results across the vertical plane we call getQfull:

fem = fe.getQfull(fem);

which computes the observation operators necessary for the interpolation of the computed quantities at the triangular elements midpoints. This provides the following additional operators, which are also stored under fem.Q:

field name	dimension	description
Qefull	nobs x ndofs	interpolates $E_x$ across the computational domain
Qhfull	nobs x ndofs	interpolates $H_x$ across the computational domain
QJefull	nobs x ndofs	interpolates $J_x$ across the computational domain
QeHyfull	nobs x ndofs	interpolates $H_y$ across the computational domain
QeHzfull	nobs x ndofs	interpolates $H_z$ across the computational domain
QhEyfull	nobs x ndofs	interpolates $E_y$ across the computational domain
QhEzfull	nobs x ndofs	interpolates $E_z$ across the computational domain

Boundary conditions

Boundary conditions are required to ensure a unique solution to the PDE. We assume a 1-D layered at both the left and right boundaries. Note that the boundaries may be different!

After successful setup of the FE discretization using FEMproblem and the assembly of the linear system using FEMassemble, the Dirichlet boundary values for all four domain boundaries are calculated in the function removeDirichlet.

A typical call would look like

fem = fe.removeDirichlet(fem, 'sigmaBCL', [1e-2, 1e-1, 1e-3], 'thicknessBCL', [2e4, 2e3], ...
        'sigmaBCL', [1e-1, 1e-3], 'thicknessBCL', 2e4, 'frequency', freq);

In the given example, the conductivities and thicknesses at the left and right boundaries are different.

At the top and bottom model boundaries, the Dirichlet values are obtained by horizontal interpolation considering a Neumann-type condition imposed at the four corners of the modelling domain.

If a simple 1-D background with identical layers at the left and the right model boundaries is considered, the two keyword pairs sigmaBCL and sigmaBCR as well as thicknessBCL and thicknessBCR can be replaced by the keyword pair sigmaBC and thicknessBC in the argument list of removeDirichlet:

fem = fe.removeDirichlet(fem, 'sigmaBC', [1e-2, 1e-3], 'thicknessBC', 2e4, freq);

Although arbitrary topography can be modelled in the computational domain, the Earth's surface at both boundaries is restricted to be at $z=0$ m. This feature will be included in a future release of the FEMT2D.

Solution of the system

The assembled system of linear equations can be solved using the call

sol = fe.FEMsolve(fem, 'verbose', true);

Depending on the desired polarization, the MATLAB structure sol contains the degrees of freedom (DOFS) sol.ue for E-polarization and/or sol.uh for the H-polarization.

Note that the actual values at the desired points of interest have to be interpolated using the afformentioned Q operators.

To obtain all possible field components for both polarizations at the positions defined by obs and one frequency fem.app.frequency, one would call

iw = -2i * pi * fem.app.frequency
Ex = fem.Q.Qe   * sol.ue;
Ey = fem.Q.QhEy * sol.uh;
Ez = fem.Q.QhEz * sol.uh;
Hx = fem.Q.Qh   * sol.uh;
Hy = fem.Q.QeHy * sol.ue / iw;
Hz = fem.Q.QeHz * sol.ue / iw;

The last two statements originate from the componentwise evaluation of Maxwell's curl equation $$\nabla \times \mathbf E = -i\omega\mu \mathbf H.$$ Note that FEMT2D is able to incorporate inhomogenous magnetic permeabilities, which are taken into account in the forming of fem.Q.QeHy and fem.Q.QeHz when calling getQ.

Visualization of the results

A plot of the modulus of the electrical current density can be produced using a call to patchplotConst:

plot.patchplotConst(fem.mesh.node, ...
    [fem.elem2dofs(1:3, :); tools.asRow(abs(fem.Q.QJefull * sol.ue))], ...
    @(x) x, 'none');
colormap(gca, jet(256));
h = colorbar();
ylabel(h, '|J| in A/m^2');
axis equal tight ij

patchplotConst takes in its first argument a list of nodes as obtained by triangle. The second argument is an 3 x nt FE element table augmented with an array of scalar quantities in an additional 4th row. In our example, this scalar quantitiy is the absolute value of the induced current density, $|J_x|$, for E-polarization. To transform this quantity, e.g., by using its common logarithm, we apply a function defined by its function handle.

Examples of function handle	Description
@(x) x	no transformation is applied, i.e., the data will be displayed as is
@log10	the common logarithm is applied to the data

The last argument controls the color of the triangular mesh. Colors can be defined either as normalized RGB values, e.g., [0.1, 0.1, 0.1] or as color abbreviation, e.g., 'w' for white. The mesh can be suppressed by using 'none' as mesh color.

The triangular elements are colored according to the color table set by a call to, e.g., colormap(gca, jet(256)).

All plot routines have been called by the top level function plot.plotMT() which accepts the following parameters:

Parameter	Value	Purpose
fem	fem struct
sol	sol struct with results of FEMsolve
mesh	mesh struct
obs	2 x np array with np observation points
periods	nf array of periods	required for sounding curve plots
profile	rhoa+phase, tipper	plot for a single frequency along profile defined by obs
sounding	rhoa+phase	sounding curves at a particular station given by its number site
section	Jx	depth section of a given field component, here $|J_x|$
section	conductivity	depth section of conductivities
section	subdomain	depth section of subdomain number which correspond to entries of the sigma and mu arrays required in the call to FEMproblem

More plot types will be added to +plot/plotMT.m in the future.

License

MIT License

Copyright (C) 2020 Ralph-Uwe Börner

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

femt2d's Issues

Octave compatible?

Compatible to free matlab clone octave?
No license is needed for octave.

If yes, please extend your info.
if unknown, untested with octave is also nice.
if no, it is also ok. Some commands are not compatible in octave.

So you are the expert. Perhaps you can make a fast test run.
Here some last infos about octave.
https://wiki.octave.org/Differences_between_Octave_and_Matlab
https://wiki.octave.org/Release_History
https://octave.org/NEWS-7.html

How to creat *.node ,*.ele,and *.1.ploy file?

Sir,I modified your notebook code to create my mesh of model,but I do not know how to creat *.node , *.ele, and *.1.ploy files,Please help me,thank a lot!

mt.exportPLC(geom,
'../meshes/layered.poly',
float_format='.8e')

print(geom)

callTriangle('../meshes/layered.poly',
quality=34.2,
verbose=False)

Is there a fast way to calculate a Jacobian (sensitivity) matrix in MT 2-D inversion with your code?

Excuse me, sir, i have a doubt?

There is a difference between the analytical solution and driverlayered.m about the layered model
This picture is layer model

This Python code is analytical method

import math
import cmath
import time


print('====================================');
print('1D MAGNETOTELLURIC MODELLING PROGRAM');
print('====================================');
print('    LAST UPDATED 17TH DECEMBER 2013 ');
print('    DEVELOPED BY ANDREW PETHICK     ');
print('      WWW.DIGITIALEARTHLAB.COM      ');
print('====================================');
print('');
print('licensed under WTFPL')
print('');
start = time.clock();

mu = 4*math.pi*1E-7; #Magnetic Permeability (H/m)
resistivities = [100, 1100, 600,10];
thicknesses = [1000,2000,3400];
n = len(resistivities);

import numpy as np
l = np.logspace(-2,3,24)
sl = (1/l).tolist()

frequencies = sl;

print('freq     \t ares   \t\t\t   phase');
for frequency in frequencies:   
    w =  2*math.pi*frequency;       
    impedances = list(range(n));
    #compute basement impedance
    impedances[n-1] = cmath.sqrt(w*mu*resistivities[n-1]*1j);


   
    for j in range(n-2,-1,-1):
        resistivity = resistivities[j];
        thickness = thicknesses[j];
  
        # 3. Compute apparent resistivity from top layer impedance
        #Step 2. Iterate from bottom layer to top(not the basement) 
        # Step 2.1 Calculate the intrinsic impedance of current layer
        dj = cmath.sqrt((w * mu * (1.0/resistivity))*1j);
        wj = dj * resistivity;
        # Step 2.2 Calculate Exponential factor from intrinsic impedance
        ej = cmath.exp(-2*thickness*dj);                     
    
        # Step 2.3 Calculate reflection coeficient using current layer
        #          intrinsic impedance and the below layer impedance
        belowImpedance = impedances[j + 1];
        rj = (wj - belowImpedance)/(wj + belowImpedance);
        re = rj*ej; 
        Zj = wj * ((1 - re)/(1 + re));
        impedances[j] = Zj;    

    # Step 3. Compute apparent resistivity from top layer impedance
    Z = impedances[0];
    absZ = abs(Z);
    apparentResistivity = (absZ * absZ)/(mu * w);
    phase = math.atan2(Z.imag, Z.real)
    print(frequency, '\t', apparentResistivity, '\t', phase);
    
print('');
print('time taken = ', time.clock() - start, 's');

I modified your driverlayered.m to satisfy this model

clearvars();
close('all');
clc();

%%
meshFile = 'meshes/layeredz1.1';

mesh = mesh.getMesh('filename', meshFile, ...
    'format', 'triangle', ...
    'shift', 0, 'scale', 1, 'verbose', true);

%%
figure(1);
plot.plotMT('mesh', mesh, 'section', 'subdomains');

%%
periods = logspace(-2, 3, 24);
freq = 1 ./ periods;
omega = 2 * pi * freq;
mu0 = pi * 4e-7;

sigma = 1.0 ./ [10, 600, 1100, 100, 1e9];
mu = [1, 1, 1, 1,1];
sigmaBCL = 1.0 ./ [100,1100, 600,10];
muL = [1, 1, 1, 1];
thkBCL = [1e3, 2e3,3.4e3];

xobs = tools.asRow(-1e4:1e3:1e4);
obs = [xobs; zeros(size(xobs)) + 0.1];

%%
fem = fe.FEMproblem('mesh', mesh, ...
    'elementtype', 'Lagrange', ...
    'order', 2, 'dimension', 2, ...
    'sigma', sigma, 'mu', mu, ...
    'application', 'MT', ...
    'polarization', 'both', ...
    'frequency', freq(1), ...
    'verbose', true);

%%
figure(2);
plot.plotMT('fem', fem, 'section', 'conductivity');

%%
fem = fe.getQ(fem, obs);
% fem = fe.getQfull(fem);

fem = fe.FEMassemble(fem, 'output', 'matrices', ...
    'verbose', false);

solution.rhoaxy = zeros(length(xobs), length(periods));
solution.rhoayx = zeros(length(xobs), length(periods));
solution.phixy = zeros(length(xobs), length(periods));
solution.phiyx = zeros(length(xobs), length(periods));


for k = 1:length(freq)
    freq = 1.0 / periods(k);

    
    fem = fe.removeDirichlet(fem, ...
        'sigmaBC', sigmaBCL, ...
        'thicknessBC', thkBCL, ...
        'frequency', freq);

    sol = fe.FEMsolve(fem, 'verbose', false);

    sol = mt.postProcessing(fem, sol);
    solution.rhoaxy(:, k) = sol.rhoaxy;
    solution.rhoayx(:, k) = sol.rhoayx;
    solution.phixy(:, k) = sol.phixy;
    solution.phiyx(:, k) = sol.phiyx;
end

%%
figure(3)
plot.plotMT('fem', fem, 'sol', solution, ...
    'sounding', 'rhoa+phase', ...
    'station', 12, 'periods', periods);

I modified your Layered_Halfspace.ipynb to mesh this model

import pygimli as pg
import matplotlib.pyplot as plt
import pygimli.meshtools as mt
import os
from os import system
import numpy as np

                 quality=33, 
                 triangle='triangle',
                 verbose=True):
    
    filebody = filename.replace('.poly', '')

    syscal = triangle + ' -pq' + str(quality)
    syscal += 'Aa ' + filebody + '.poly'
 
    if verbose:
        print(syscal)

    system(syscal)

world = mt.createWorld(start=[-4e5, -4e5], 
                       end=[4e5, 4e5], 
                       layers=[0, 1e3, 3e3,6.4e3], 
                       area=[0, 1e7, 1e6, 1e6,0],
                       marker=[1, 2, 3, 4,5],
                       worldMarker=False)

geom = world

ax, _ = pg.show(geom, 
                showNodes=False, 
                boundaryMarker=False)
ax.set_ylim(ax.get_ylim()[::-1]);

mt.exportPLC(geom,
             '../meshes/layeredz1.poly', 
             float_format='.8e')

print(geom)

callTriangle('../meshes/layeredz1.poly', 
             quality=34.2,
             verbose=False)

mesh = mt.createMesh(geom, quality=34.2)
ax, _ = pg.show(mesh)
ax.set_ylim(ax.get_ylim()[::-1]);

Then I use ./triangle -pq34.2Aa layeredz1.poly in the terminal to generate layeredz1.1.ele,layeredz1.1.node,layeredz1.1.poly.
this zip file contain these four files.
layeredz1.zip

this is analytical solution

this solution is after modifying the driverlayered.m