ruboerner / femt2d Goto Github PK

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');

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);

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]);