Influence of floating point precision on drift-diffusion simulation

This note­books shows the nu­mer­ic ef­fects in fi­nite vol­ume drift-dif­fuion cal­cu­la­tions. A par­tic­u­lar­ly dif­fi­cult case of com­pu­ta­tion is cho­sen, where charge car­ri­er den­si­ty is very high in­side the de­vice, but the den­si­ty of cur­rent is low. More­o­ev­er, com­pu­ta­tion­al grid is high­ly re­fined close to the elec­trode to fur­ther ex­ager­ate the nu­mer­ic ef­fects.

The sam­ple cal­cu­la­tion is run us­ing with var­ied pre­ci­sion, and the re­sults are com­pared. The choice of pre­ci­sion are:

  • dou­ble pre­ci­sion (np.double), with ~15 sig­nif­i­cant dig­its
  • ex­tend­ed pre­ci­sion (np.float128), with ~18 sig­nif­i­cant dig­its
  • if pos­si­ble, with quadru­ple pre­ci­sion (oedesext.precision.qfloat), with ~30 sig­nif­i­cant dig­its

Note that numpy cur­rent­ly np.float128 us­es on­ly 80 bits to store the val­ues, and the re­main­ing bits are not used.

import numpy as np
types = [np.double, np.float128]
try:
    # oedesext are extensions to oedes not yet released as open-source
    import oedesext.precision
    types.append(oedesext.precision.qfloat)
except:
    oedesext=None
%matplotlib inline
from oedes import *
import matplotlib.pylab as plt
import scipy.constants
# Set up a hole only device with built-in voltage of 1V, high charge carrier
# density on one side, minimum mesh spacing of 0.1 Å, biased with 1 mV.
b = models.BaseModel()
models.electronic_device(b, fvm.mesh1d(100e-9, dx_boundary=1e-11), 'p')
b.poisson.bc = [ models.AppliedVoltage('electrode0',calculate_current=True),
                 models.AppliedVoltage('electrode1',calculate_current=True) ]
b.setUp()
b.species[0].convergenceTest = fvm.ElementwiseConvergenceTest(
    rtol=1e-12, atol=0)
params = {
    'T': 300,
    'electrode0.voltage': 1e-3,
    'electrode0.workfunction': 0,
    'electrode1.voltage': 0,
    'electrode1.workfunction': -1,
    'hole.mu': 1e-9,
    'hole.N0': 1e27,
    'hole.energy': 0,
    'epsilon_r': 3}

Residuals

The norm of resid­u­als rough­ly cor­re­sponds to the float­ing point pre­ci­sion.

sol = {}
for t in types:
    sol[t] = solve(b, np.asarray(b.X, dtype=t), params, maxiter=30)
    residuals = b.residuals(0., sol[t], 0. * sol[t],params)
    print('%20s |F|=%e' %(t,np.linalg.norm(residuals)))

Numeric estimation of charge carrier density

All so­lu­tions are the same.

for t in types:
    o = b.output(0., sol[t], 0. * sol[t], params)
    plt.plot(b.poisson.mesh.cells['center'],np.abs(o['hole.c']),label=t)
plt.yscale('log')
plt.legend(loc=0);

Numeric current density

De­spite that all so­lu­tions are cor­rect, not all es­ti­ma­tions of cur­rent are cor­rect:

  • es­ti­ma­tion electrode1.Jboundary, based on cur­rent den­si­ty at x=100 nm, is cor­rect for all pre­ci­sions
  • es­ti­ma­tion electrode0.Jbounary, based on cur­rent den­si­ty at x=0 nm, is wrong for all ex­cept pre­ci­sions ex­cept for quadru­ple pre­ci­sion
  • av­er­aged es­ti­ma­tion J, based on Ramo-Shock­ley the­o­rem, is mod­er­ate­ly ac­cu­rate, and, in this ex­am­ple, suf­fi­cient in ex­tend­ed pre­ci­sion
for t in types:
    out = b.output(0, sol[t], 0. * sol[t], params)
    d=dict(t=t,J=out['J'],Jb0=out['electrode0.Jboundary'],
           Jb1=out['electrode1.Jboundary'])
    print('{t!s} J={J:e} J_boundary={Jb0:e},{Jb1:e}'.format(**d))

Numeric current continuity

Be­low nu­mer­ic cur­rent con­ti­nu­ity close to the left elec­trode is plot­ted for dif­fer­ent choic­es of float­ing point pre­ci­sion. Sur­pris­ing­ly, good charge car­ri­er den­si­ty can be cal­cu­lat­ed de­spite of very es­ti­ma­tion of cur­rent den­si­ty.

for t in types:
    o = b.output(0., sol[t], 0. * sol[t], params)
    j = o['hole.j'] * b.species[0].z * scipy.constants.elementary_charge
    plt.plot(b.poisson.mesh.faces['center'],np.abs(j),label=t)
plt.yscale('log')
plt.xscale('log')
plt.xlabel('x [m]')
plt.ylabel('j [A/m2]')
plt.legend(loc=0);