<jupyter_start><jupyter_text>DC Resistivity Forward SimulationThis example runs a [direct-current resistivity sounding](https://gpg.geosci.xyz/content/DC_resistivity/index.html) over a layered earth. The array configuration is a [schlumberger array](https://gpg.geosci.xyz/content/DC_resistivity/DC_surveys.html?highlight=schlumberger). The potential electrodes are centered between the closest-spaced current electrodes. If you run into any problems, please report an [issue](https://github.com/lheagy/gwb-dc-dashboard/issues/new). The source code is available on github [DC-Forward-Simulation.ipynb](https://github.com/lheagy/gwb-dc-dashboard/blob/master/DC-Forward-Simulation.ipynb)<jupyter_code># %matplotlib widget
# import ipympl

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import discretize 
import ipywidgets as widgets
from SimPEG import (
    DC, Maps, DataMisfit, Regularization,
    Optimization, Inversion, InvProblem, Directives, Utils
)
from SimPEG.EM.Static.Utils.StaticUtils import apparent_resistivity
from pymatsolver import Pardiso

# plt.ioff()
from matplotlib import rcParams
rcParams['font.size'] = 12
# survey parameters

def schlumberger_array(min_spacing=2., max_spacing=250., n_soundings=20, z_loc=-0.01):
    
    a0 = np.log10(min_spacing/2.)
    an = np.log10(max_spacing/2.)
    
    loc_a = np.logspace(a0, an, n_soundings)
    loc_b = loc_a.copy()
    loc_a = -loc_a

    loc_m = loc_a[0]/2.
    loc_n = loc_b[0]/2.

    z_src = -0.01
    
    return (
        np.vstack([loc_a, np.ones_like(loc_a)*z_src]).T,
        np.vstack([loc_b, np.ones_like(loc_b)*z_src]).T,
        np.atleast_2d(np.hstack([loc_m, z_src])),
        np.atleast_2d(np.hstack([loc_n, z_src])),
    )
# create a mesh

def simulation_mesh(electrode_spacing, n_cells_per_electrode=4, padding_factor=1.3, padding_distance_factor=1.5): 
    dx = np.diff(electrode_spacing) if len(electrode_spacing) > 1 else electrode_spacing
    core_h = np.kron(dx/n_cells_per_electrode, np.ones(n_cells_per_electrode))

    def add_padding(hx, pad_factor=padding_factor, pad_distance=np.max(electrode_spacing)*padding_distance_factor):
        while(sum(hx)) < pad_distance:
            pad = hx[-1]*pad_factor
            hx = np.hstack([hx, np.r_[pad]])
        return hx

    h = add_padding(core_h)
    hz = np.flipud(h.copy())

    hx = np.hstack([np.flipud(h), h])
    mesh = discretize.TensorMesh([hx, hz], x0='CN')
    
    return mesh
def inversion_mesh(simulation_mesh):
    return discretize.TensorMesh([simulation_mesh.hy], x0='N')

def synthetic_model(inversion_mesh, layer_tops, layer_rho):
    model = np.zeros(inversion_mesh.nC)

    for t, r in zip(layer_tops, layer_rho):
        model[inversion_mesh.gridCC <= t] = r
    
    return model    
# set up the survey and run a forward simulation

def create_simulation(mesh, mapping, loc_a, loc_b, loc_m, loc_n):

    src_list = [
        DC.Src.Dipole(
            [DC.Rx.Dipole_ky(loc_m, loc_n)], 
            loc_a[i, :], loc_b[i, :]
        ) for i in range(loc_a.shape[0])
    ]

    survey = DC.Survey_ky(src_list)
    prob = DC.Problem2D_N(mesh, rhoMap=mapping, storeJ=True, Solver=Pardiso)
    prob.pair(survey)
    
    return prob, survey
def staircase_model(x, mod): 
    new_x = np.hstack([np.kron(x[:-1], np.r_[1, 1])[1:], np.r_[x[-1]]])
    new_m = np.kron(mod, np.r_[1, 1]) 
    
    return new_x, new_m
def synthetic_data(layer_tops, layer_rho, min_spacing=2., max_spacing=250., n_soundings=20):

    loc_a, loc_b, loc_m, loc_n = schlumberger_array(min_spacing, max_spacing, n_soundings)
    electrode_spacing = loc_b[:, 0] - loc_a[:, 0]
    mesh = simulation_mesh(electrode_spacing)
    inv_mesh = inversion_mesh(mesh)
    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh)
    true_model = np.log(synthetic_model(inv_mesh, layer_tops, layer_rho))

    prob, survey = create_simulation(mesh, mapping, loc_a, loc_b, loc_m, loc_n)

    fields = prob.fields(true_model)
    dobs = survey.dpred(true_model, f=fields)

    survey.dobs = dobs

    return prob, survey, fields
def plot_model_data(
    prob, survey, fields=None, view="model", source_ind=None, ax=None, xlim=None, zlim=None
):
    
    if ax is None: 
        fig, ax = plt.subplots(1, 2, figsize =(14, 5))
    
    locs = np.vstack([np.r_[src.loc[0][0], src.loc[1][0]] for src in survey.srcList])
    electrode_spacing = locs[:, 1] - locs[:, 0]
    max_electrode_spacing = np.max(locs)
    
    if xlim is None: 
        xlim = max_electrode_spacing * 1.1 * np.r_[-1, 1]
    if zlim is None: 
        zlim = np.r_[-max_electrode_spacing, np.min(np.abs(locs))*4]
    
    if view.lower() == "model": 
        cb = plt.colorbar(
            prob.mesh.plotImage(
                prob.rho, ax=ax[0], pcolorOpts={'norm':LogNorm()}, showIt=False
            )[0], 
            ax=ax[0]
        )
        cb.set_label("Resistivity ($\Omega$m)")
        ax[0].set_title("model")
        
    else:
        src = survey.srcList[source_ind]
        phi = DC.Rx.IntTrapezoidal(prob.kys, fields[src, "phi"])
        if view.lower() == "potential":
            phi = prob.mesh.aveN2CC * phi
            cb = plt.colorbar(
                prob.mesh.plotImage(phi, ax=ax[0], showIt=False)[0], 
                ax=ax[0]
            )
            cb.set_label("Electric potential (V/m)")
            ax[0].set_title("potential")
        elif view.lower() == "currents":
            j = prob.mesh.aveE2CCV * (prob.MeI * (prob.MeSigma * (- prob.mesh.nodalGrad * phi)))
            cb = plt.colorbar(
                prob.mesh.plotImage(
                    j, ax=ax[0], showIt=False, view='vec', vType='CCv', pcolorOpts={'norm':LogNorm()},
                    streamOpts={'color':'k', 'linewidth': 1}, 
                    stream_threshold=1e-5*np.max(np.abs(j)),
                    range_x=xlim, range_y=zlim, sample_grid=[(xlim[1]-xlim[0])/50, (zlim[1]-zlim[0])/50]
                )[0], 
                ax=ax[0]
            )
            cb.set_label("Current density (A/m$^2$)")
            ax[0].set_title("current density")
        
        
    max_wire_height = zlim[1]*0.75
    wire_height = max_wire_height/locs.shape[0]
    
    if source_ind is not None:
        locs = np.atleast_2d(locs[source_ind, :])
        
    for i in range(locs.shape[0]):  
        wirex = np.kron(locs[i, :], np.r_[1, 1])
        wirez = np.r_[0, wire_height*(i+1), wire_height*(i+1), 0]
        ax[0].plot(wirex, wirez, "k", lw=0.5, alpha=0.5)

    ax[0].plot(locs[:,0], np.zeros_like(locs[:,0]), 'C1v', ms=4)
    ax[0].plot(locs[:,1], np.zeros_like(locs[:,1]), 'C1v', ms=4)
    ax[0].set_aspect(1)

    ax[0].set_xlim(xlim)
    ax[0].set_ylim(zlim)
    
    rho_a = apparent_resistivity(survey, dobs=survey.dobs)
    ax[1].loglog(electrode_spacing/2, rho_a, 'o')
    ax[1].set_xlabel("electrode separation / 2 (m)")
    ax[1].set_ylabel("apparent resistivity ($\Omega$m)")
    ax[1].set_title("simulated data")
    ax[1].grid(which="both", alpha=0.5, lw=0.5)
    
    if source_ind is not None:
        ax[1].loglog(electrode_spacing[source_ind]/2, rho_a[source_ind], 'C1o', ms=8)
    
    rholim = [rho_a.min()*0.5, rho_a.max()*2]
    ax[1].set_ylim(rholim)

    plt.tight_layout()
    
    return ax
class ForwardSimulationApp:
    
    def __init__(self):
        
        # Model Parameters
        self.nlayer = widgets.IntSlider(min=1, max=5, description="N layers")
        self.layertops = widgets.HBox([self._layer_top_widget(i) for i in range(self.nlayer.value)])
        self.layerrho = widgets.HBox([self._layer_rho_widget(i) for i in range(self.nlayer.value)])
        
        self.layertops.children[0].disabled = True

        self.nlayer.observe(self._update_box_layout, names=["value"])

        model_title = widgets.HTML(
            value="""
            <h2>Model Parameters</h2>
            <p>
            These describe the setup of the earth model. You can add up to 5 layers. 
            For each layer, define the depth to the top of the layer (in meters) and 
            the resistivity (in &Omega; m)
            </p>
            """,
        )

        w_model = widgets.VBox([model_title, self.nlayer, self.layertops, self.layerrho])
        
        # Survey Parameters
        self.min_spacing = widgets.BoundedFloatText(value=2., min=0., description="min $AB$")
        self.max_spacing = widgets.BoundedFloatText(value=100., min=0, max=1000, description="max $AB$")
        self.n_soundings = widgets.IntSlider(min=1, max=50, value=20, description="N soundings")
        survey_title = widgets.HTML(
            value="""
            <h2>Survey Parameters</h2>
            <p>
            Define the source electrode locations. The receiver electrodes are 
            positioned at half of the distance between the smallest source electrode 
            spacing. For example, if the minimum source electrode spacing is 2m, then 
            the receiver electrodes (M and N electrodes) dipole will be from -0.5m and 0.5m. 
            </p>
            """,
        )
        
        self.min_spacing.observe(self._update_limits, names=["value"])
        self.max_spacing.observe(self._update_limits, names=["value"])
        self.n_soundings.observe(self._recompute, names=["value"])
    
        w_survey = widgets.VBox([survey_title, self.n_soundings, self.min_spacing, self.max_spacing])
    
        # View parameters
        self.view_toggles = widgets.ToggleButtons(
            options=["model", "potential", "currents"],
            tooltips=["resistivity of the earth", "electric potential"],
            value="model"
        )
        self.view_toggles.observe(self._update_view, names=["value"])
        self.view = widgets.HBox([self.view_toggles])
        self._xlim_slider()
        self._zlim_slider()
        limits = widgets.Box([widgets.VBox([self.xlim, self.zlim])]) #, self.rholim])
            
        self.source_slider = None
        view_title = widgets.HTML(
            value="""
            <h2>View Parameters</h2>
            """
        )
        
        w_view = widgets.VBox([view_title, self.view, limits])
        
        # Assemble widgets
        w = widgets.VBox([w_model, w_survey, w_view])
        
        # plot container
        self._plot_container = widgets.Output()
        self.container = widgets.VBox([w, self._plot_container])
        
        self._update_app()
    
    def __call__(self):
        return self.container
    
    def _layer_top_widget(self, i):
        w = widgets.BoundedFloatText(value=10*i, min=0, max=1000, description=f"depth_{i}", layout=widgets.Layout(width="20%"))
        w.observe(self._update_limits, names=["value"])
        return w
        
    def _layer_rho_widget(self, i):
        w = widgets.BoundedFloatText(value=100, min=1e-8, max=1e8, description="{}_{}".format("$\\rho$", i), layout=widgets.Layout(width="20%"))
        w.observe(self._update_limits, names=["value"])
        return w
    
    def _update_box_layout(self, change):
        new = change["new"]
        old = change["old"]
        if new < old:
            self.layertops.children = self.layertops.children[:new]
            self.layerrho.children = self.layerrho.children[:new]
        elif new > old:
            self.layertops.children += tuple([self._layer_top_widget(i+old) for i in range(new-old)])
            self.layerrho.children += tuple([self._layer_rho_widget(old+i)for i in range(new-old)])
        self._update_limits(change)
    
    def _update_view(self, change):
        new = change["new"]
        old = change["old"]
        if new == "model":
            self.view.children = self.view.children[:1]
            self.source_slider = None
        else:
            self.source_slider = widgets.IntSlider(
                value=0, min=0, max=self.n_soundings.value-1, description="source #"
            )
            self.source_slider.observe(self._redraw, names=["value"])
            if len(self.view.children) == 1:
                self.view.children += (self.source_slider,)
        self._redraw(None)

    def _xlim_slider(self):
        max_electrode_spacing = self.max_spacing.value/2
        
        if getattr(self, 'xlim', None) is None:
            xlim = widgets.FloatRangeSlider(
                value=tuple(max_electrode_spacing * 1.1 * np.r_[-1, 1]),
                min=-max_electrode_spacing*2,
                max=max_electrode_spacing*2,
                description="xlim",
                continuous_update=False,
                readout_format='.0f'
            )
            self.xlim = xlim
            self.xlim.observe(self._redraw, names="value")
        else:
            self.xlim.min = -max_electrode_spacing*2
            self.xlim.max = max_electrode_spacing*2
    
    def _zlim_slider(self):
        max_electrode_spacing = self.max_spacing.value/2
        min_electrode_spacing = self.min_spacing.value/2
        if getattr(self, 'zlim', None) is None:
            zlim = widgets.FloatRangeSlider(
                value=[-max_electrode_spacing*1.5, min_electrode_spacing*4],
                min=-max_electrode_spacing*2,
                max=min_electrode_spacing*8,
                description="zlim",
                continuous_update=False,
                readout_format='.0f'
            )
            self.zlim = zlim
            self.zlim.observe(self._redraw, names="values")
        else:
            self.zlim.min = -max_electrode_spacing*2
            self.zlim.max = min_electrode_spacing*8
            
    
    def _update_limits(self, _):
#         self._rholim_slider()
        self._xlim_slider()
        self._zlim_slider()
        self._update_app(True)

    def _recompute(self, _):
        self._update_app(recompute=True)
        
    def _redraw(self, _):
        self._update_app(recompute=False)
        
    def _create_plot(self, recompute=True):
        fig, ax = plt.subplots(1, 2, figsize =(11, 3.5), dpi=120)
        layer_tops = -1*np.array([top.value for top in self.layertops.children])
        layer_rho = np.array([rho.value for rho in self.layerrho.children])
        
        if recompute is True or getattr(self, '_fields', None) is None:
            prob, survey, fields = synthetic_data(
                layer_tops, layer_rho, self.min_spacing.value, self.max_spacing.value, self.n_soundings.value
            )
            self._prob = prob
            self._survey = survey
            self._fields = fields
            
        plot_model_data(
            self._prob, self._survey, self._fields, view=self.view_toggles.value, 
            source_ind=None if self.source_slider is None else self.source_slider.value,
            ax=ax, xlim=self.xlim.value, zlim=self.zlim.value
        )
        
    def _update_app(self, recompute=True):
        self._plot_container.clear_output(wait=True)
        with self._plot_container:
            self._create_plot(recompute)
            plt.show()
        
app = ForwardSimulationApp()
app()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12956942877993763, 'fraction_numerical': 0.015325416307304452, 'mean_word_length': 2.4549621819848726, 'pattern_counts': {'":': 4, '<': 18, '<?xml version=': 0, '>': 17, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '522518', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5375, 'n_tokens_neox': 5143, 'n_words': 938, 'path': 'DC-Forward-Simulation.ipynb', 'repo_name': 'lheagy/gwb-dc-dashboard', 'size': '15073'}