Mesh: Basic Forward 2D DC ResistivityΒΆ

2D DC forward modeling example with Tensor and Curvilinear Meshes

../../../_images/sphx_glr_plot_poisson_001.png
from SimPEG import Mesh, Utils, SolverLU
import numpy as np
import matplotlib.pyplot as plt


def run(plotIt=True):
    # Step1: Generate Tensor and Curvilinear Mesh
    sz = [40, 40]
    tM = Mesh.TensorMesh(sz)
    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz, 'rotate'))

    # Step2: Direct Current (DC) operator
    def DCfun(mesh, pts):
        D = mesh.faceDiv
        sigma = 1e-2*np.ones(mesh.nC)
        MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
        A = -D*MsigI*D.T
        A[-1, -1] /= mesh.vol[-1]  # Remove null space
        rhs = np.zeros(mesh.nC)
        txind = Utils.meshutils.closestPoints(mesh, pts)
        rhs[txind] = np.r_[1, -1]
        return A, rhs

    pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5]))

    # Step3: Solve DC problem (LU solver)
    AtM, rhstM = DCfun(tM, pts)
    AinvtM = SolverLU(AtM)
    phitM = AinvtM*rhstM

    ArM, rhsrM = DCfun(rM, pts)
    AinvrM = SolverLU(ArM)
    phirM = AinvrM*rhsrM

    if not plotIt:
        return

    # Step4: Making Figure
    fig, axes = plt.subplots(1, 2, figsize=(12*1.2, 4*1.2))
    vmin, vmax = phitM.min(), phitM.max()
    dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
    dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
    cb0 = plt.colorbar(dat[0], ax=axes[0])
    cb1 = plt.colorbar(dat[0], ax=axes[1])
    cb0.set_label("Voltage (V)")
    cb1.set_label("Voltage (V)")

    axes[0].set_title('TensorMesh')
    axes[1].set_title('CurvilinearMesh')


if __name__ == '__main__':
    run()
    plt.show()

Total running time of the script: ( 0 minutes 2.187 seconds)

Generated by Sphinx-Gallery