Attempt to solve MOSCAP problem posed in https://github.com/usnistgov/fipy/issues/731

MOSCAP.ipynb

I realise this is a very old thread, but in case anyone is interested in testing this further, I believe this semiconductor system can be solved quite decently using Gummel iterations, similar to the strategy used in https://github.com/nsdt-zhaw/ChargeFabrica/. One important thing to note is that pixel width should be at or below 1 nm for accurate results. Also, because electron density for some pixels goes above 1e25, I think some entries in the residual vector completely dominate, which makes it hard to tell exactly how many iterations are needed for "good" convergence, but from my testing it looks like the density values are physical (above zero) up to 1V (I haven't tested above).

Simulation Code File:

# -*- coding: utf-8 -*-
#This code is a MOSCAP simulation
import os
os.environ["OMP_NUM_THREADS"] = "1" #Really important! Pysparse doesnt benefit from multithreading.
import numpy as np
from fipy import CellVariable, TransientTerm, DiffusionTerm, ExponentialConvectionTerm, ImplicitSourceTerm, ResidualTerm
import fipy
from fipy.tools import numerix
import time
from scipy import constants
from joblib import Parallel, delayed
import multiprocessing
import copy

TInfinite = 300.0
(q, epsilon_0, D) = (constants.electron_volt, constants.epsilon_0, (constants.Boltzmann * TInfinite) / constants.electron_volt)

def solve_for_voltage(voltage):
    solver = fipy.solvers.pysparse.linearLUSolver.LinearLUSolver(precon=None, iterations=1, tolerance=1.00e-10) #Very fast solver

    L = 5e-6
    h = 10e-6
    tox = 0.1 * 1e-6
    q = 1.6 * 1e-19
    un = 0.14
    up = 0.045
    p0 = 1e21
    n0 = 1e11
    e0 = 8.854 * 1e-12
    nx = 1

    mesh1 = fipy.Grid1D(dx=h / 10000, nx=10000)  # mesh for domain 1 for solving for p and n
    mesh2 = fipy.Grid1D(dx=tox / 10, nx=10)
    mesh3 = mesh1

    y = mesh3.cellCenters[0]

    N = 1e21 * (y <= h)  # N changes from Domain 1 to Domain 2
    e = 11.9 * e0 * (y <= h) + 3.9 * e0 * (y > h)  # e changes from Domain 1 to Domain 2

    philocal = CellVariable(name="electrostatic potential", mesh=mesh3, value=0.00, hasOld=True)
    nlocal = CellVariable(name="electron density", mesh=mesh3, value=n0, hasOld=True)
    plocal = CellVariable(name="hole density", mesh=mesh3, value=p0, hasOld=True)

    philocal.constrain(0., where=mesh3.facesLeft)
    philocal.constrain(voltage, where=mesh3.facesRight)
    plocal.constrain(p0, where=mesh3.facesLeft)
    plocal.constrain(p0 * numerix.exp(-philocal.faceValue / D), where=mesh3.facesRight)
    nlocal.constrain(n0, where=mesh3.facesLeft)
    nlocal.constrain(n0 * numerix.exp(philocal.faceValue / D), where=mesh3.facesRight)

    eqn = (0.00 == -TransientTerm(coeff=q, var=nlocal) + DiffusionTerm(coeff=q * D * un,var=nlocal) - ExponentialConvectionTerm(coeff=q * un * philocal.faceGrad,var=nlocal) + ImplicitSourceTerm(coeff=q/1e-6,var=nlocal) - q*n0/1e-6)
    eqp = (0.00 == -TransientTerm(coeff=q, var=plocal) + DiffusionTerm(coeff=q * D * up,var=plocal) + ExponentialConvectionTerm(coeff=q * up * philocal.faceGrad,var=plocal) + ImplicitSourceTerm(coeff=q/1e-6,var=plocal) - q*p0/1e-6)
    eqpoisson = (0.00 == -TransientTerm(var=philocal) + DiffusionTerm(coeff=e, var=philocal) + q*(plocal - nlocal - N)*(y <= h) + 0*(y>h))

    dt, MaxTimeStep, desired_residual, DampingFactor, NumberofSweeps, min_iterations, max_iterations = 1e-6, 1e-5, 1e-15, 0.2, 1, 100, 2000
    residual, residual_old, dt_old, TotalTime, SweepCounter = 1., 1e10, dt, 0.0, 0

    while SweepCounter < min_iterations or (SweepCounter < max_iterations and residual > desired_residual):

        t0 = time.time()

        for i in range(NumberofSweeps):
            eqpoisson.sweep(dt = dt, solver=solver)
            philocal.setValue(DampingFactor * philocal + (1 - DampingFactor) * philocal.old) # The potential should be damped BEFORE passing to the continuity equations!

            residual = eqn.sweep(dt = dt, solver=solver) + eqp.sweep(dt = dt, solver=solver)
            nlocal.setValue(DampingFactor * np.maximum(nlocal, 1.00e-30) + (1 - DampingFactor) * nlocal.old)
            plocal.setValue(DampingFactor * np.maximum(plocal, 1.00e-30) + (1 - DampingFactor) * plocal.old)

        PercentageImprovementPerSweep = (1 - (residual / residual_old) * dt_old / dt) * 100

        if residual > residual_old * 1.2:
            dt = max(1e-9, dt * 0.1)
        else:
            dt = min(MaxTimeStep, dt * 1.05)

        dt_old = dt
        residual_old = residual

        #Update old
        for v in (nlocal, plocal, philocal): v.updateOld()

        TotalTime = TotalTime + dt

        print("Sweep: ", SweepCounter, "TotalTime: ", TotalTime, "Residual: ", residual, "Time for sweep: ", time.time() - t0, "dt: ", dt, "Percentage Improvement: ", PercentageImprovementPerSweep, "Damping: ", DampingFactor)
        SweepCounter += 1

    def reshapefunction(FipyFlattenedArray):
        return [np.reshape(arr, (10000, nx)) for arr in FipyFlattenedArray]

    E = -philocal.grad  #Vector Quantity

    Efield_matrix = np.reshape(E.globalValue, (E.shape[0], 10000, nx))

    (PotentialMatrix, NMatrix, PMatrix) = reshapefunction([philocal, nlocal, plocal])

    return {"NMatrix": NMatrix, "PMatrix": PMatrix, "PotentialMatrix": PotentialMatrix, "Efield_matrix": Efield_matrix, "n": nlocal.globalValue, "p": plocal.globalValue, "phi": philocal.globalValue, "ResidualMatrix": residual, "SweepCounterMatrix": SweepCounter}

def simulate_device(output_dir):
    applied_voltages = np.arange(0.0, 1.0, 0.05)

    chunk_size = len(applied_voltages)

    def append_to_npy(filename, new_data):
        new_data = np.expand_dims(new_data, axis=0)
        path = os.path.join(output_dir, filename)
        if os.path.isfile(path):
            fulldata = np.concatenate((np.load(path), new_data), axis=0)
        else:
            fulldata = new_data
        np.save(path, fulldata)

    # Process voltages in sequential chunks
    for start in range(0, len(applied_voltages), chunk_size):
        # Create a chunk of voltages to simulate in parallel
        chunk_voltages = applied_voltages[start:start + chunk_size]

        # Parallel computation within the chunk
        chunk_results = Parallel(n_jobs=chunk_size, backend="multiprocessing")(delayed(solve_for_voltage)(voltage) for voltage in chunk_voltages)

        #DeepCopy To avoid overwriting the results in next loop
        copied_result = [copy.deepcopy(r) for r in chunk_results]

        # Save dictionary of chunk_results as .npy files named after the key
        for result in copied_result:
            for key, value in result.items():
                append_to_npy(key + ".npy", value)

        #Save an array of all the voltages applied so far
        np.save(os.path.join(output_dir, "applied_voltages.npy"), applied_voltages[:start + len(chunk_voltages)])

    return copied_result

def main_workflow():
    current_file = os.path.basename(__file__)[0:-3]
    voltage_sweep_output_dir = "./Outputs/" + current_file + "/VoltageSweep"

    if not os.path.exists(voltage_sweep_output_dir):
        os.makedirs(voltage_sweep_output_dir)

    print("Starting standard voltage sweep...")
    results = simulate_device(output_dir=voltage_sweep_output_dir)
    print("Voltage sweep completed.")
    return results

# Fix for multiprocessing on Windows
if __name__ == '__main__':
    main_workflow()

Plotting Code File:

import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib.widgets import Slider
from scipy.interpolate import interp1d

# Load data
Simulation_folder = "./Outputs/Drift_Diffusion_1D_IV_05_01_2026_MOSCAP_Problem/VoltageSweep/"
NumberOfSuns = 1.00

PMatrix = np.load(Simulation_folder + "p.npy")[:]

print("PMatrix shape: ", PMatrix.shape)

NMatrix = np.load(Simulation_folder + "n.npy")[:]
PotentialMatrix = np.load(Simulation_folder + "phi.npy")[:]

applied_voltages = np.load(Simulation_folder + "applied_voltages.npy")[:]
print("applied_voltages ", applied_voltages)
ResidualMatrix = np.load(Simulation_folder + "ResidualMatrix.npy")[:]
SweepCounterMatrix = np.load(Simulation_folder + "SweepCounterMatrix.npy")[:]
print("SweepCounterMatrix", SweepCounterMatrix)

print "applied_voltages shape: ", applied_voltages.shape

print("PotentialMatrix shape: ", PotentialMatrix.shape)

titles = ['Potential (V)']
data_matrices = [PotentialMatrix]
num_plots = len(data_matrices)

fig, axs = plt.subplots(2, 4, figsize=(12, 10))
axs = axs.ravel()
cax_list = []
colorbars = []
norms = []

for i, matrix in enumerate(data_matrices):
    cax, = axs[i].plot(matrix[0]) # the comma unpacks the first element
    axs[i].set_title("{} at {:.3f} V".format(titles[i], applied_voltages[0]))
    cax_list.append(cax)
    #Flip the x axis
    axs[i].invert_xaxis()

fig.subplots_adjust(left=0.08, right=0.98, top=0.93, bottom=0.20, wspace=0.2, hspace=0.4)
slider_ax = fig.add_axes([0.2, 0.05, 0.6, 0.03])
slider = Slider(ax=slider_ax, label='Voltage [V]', valmin=applied_voltages.min(), valmax=applied_voltages.max(), valinit=applied_voltages[0])
interp_indices = interp1d(applied_voltages, np.arange(len(applied_voltages)), bounds_error=False, fill_value="extrapolate")

def update(val):


    voltage = slider.val
    # find closest index for the given voltage
    frame = int(np.clip(np.round(interp_indices(voltage)), 0, len(applied_voltages) - 1))

    for i, matrix in enumerate(data_matrices):
        axs[i].clear()
        axs[i].plot(matrix[frame])
        axs[i].set_title("{} at {:.3f} V".format(titles[i], applied_voltages[frame]))
        axs[i].invert_xaxis()

    axs[-3].clear()
    axs[-3].plot(NMatrix[frame][:], "r")
    axs[-3].plot(PMatrix[frame][:], "b")
    axs[-3].set_yscale("log")
    axs[-3].legend(["Electron", "Hole"])
    axs[-3].set_title("Charge Carrier Distribution (1/m3)")
    axs[-3].invert_xaxis()
    fig.canvas.draw_idle()

update(0) # Initial call to display the first frame
slider.on_changed(update)
plt.show()