Determination of the Electric Field at Contacts

[lmh91]

The issue came up from #279

Right now, we assign an electric field vector to each grid point of the grid of the electric potential.
The components of the vector is calculated by the differences in potential to the neighboring grid points in the respective dimensions. However, for 2D-contacts, e.g. a plane in x-y, a wrong value is determined for the component in z: E_z.

#279 is a hotfix for this issue such that the value which corresponds to the side lying in the semiconductor is assigned to the grid point. Which is okay for calculating the drift of charge carriers inside the semiconductor.

However, in general this is not correct.

I think instead of assigning electric field vectors to all grid points of the electric potential,
we should assign electric field vectors to each volume between the grid points of the electric potential.
Thus, the grid of the electric field would be smaller by one in each dimension
as its ticks would be the midpoints of the ticks of the grid of the electric potential.

This would also require changes in the interpolation of the electric field vector during the drift.

[lmh91]

For the electric field energy calculation we already determine the electric field for each voxel
by the 8 grid points of the corners of each voxel:

SolidStateDetectors.jl/src/Simulation/ElectricFieldEnergy.jl

Lines 22 to 68 in 3cb4d47

	function calculate_stored_energy(ep::ElectricPotential{T,3,CS}, ϵ::DielectricDistribution{T,3,CS}; consider_multiplicity::Bool = true) where {T <: SSDFloat, CS}
	cylindrical = CS == Cylindrical
	phi_2D = cylindrical && size(ep, 2) == 1
	ep3d = phi_2D ? get_2π_potential(ep, n_points_in_φ = 2) : _get_closed_potential(ep)
	grid = ep3d.grid
	W::T = 0
	for i3 in 1:size(grid, 3)-1
	for i2 in 1:size(grid, 2)-1
	for i1 in 1:size(grid, 1)-1
	w1, w2, w3 = voxel_widths(grid, i1, i2, i3)
	dV = voxel_volume(grid, i1, i2, i3, w1, w2, w3)
	_ϵ = ϵ.data[i1 + 1, i2 + 1, i3 + 1]
	ep000 = ep3d.data[i1 , i2 , i3 ]
	ep100 = ep3d.data[i1 + 1, i2 , i3 ]
	ep010 = ep3d.data[i1 , i2 + 1, i3 ]
	ep110 = ep3d.data[i1 + 1, i2 + 1, i3 ]
	ep001 = ep3d.data[i1 , i2 , i3 + 1]
	ep101 = ep3d.data[i1 + 1, i2 , i3 + 1]
	ep011 = ep3d.data[i1 , i2 + 1, i3 + 1]
	ep111 = ep3d.data[i1 + 1, i2 + 1, i3 + 1]
	efv1 = ( (ep100 - ep000) + (ep110 - ep010) + (ep101 - ep001) + (ep111 - ep011) ) / (4 * w1)
	efv2 = if cylindrical
	_w2 = (grid[2].ticks[i2 + 1] - grid[2].ticks[i2])
	if i1 == 1
	((ep110 - ep100)/(_w2*grid[1].ticks[i1+1]) +
	(ep111 - ep101)/(_w2*grid[1].ticks[i1+1])) / 2
	else
	((ep010 - ep000)/(_w2*grid[1].ticks[i1]) +
	(ep110 - ep100)/(_w2*grid[1].ticks[i1+1]) +
	(ep011 - ep001)/(_w2*grid[1].ticks[i1]) +
	(ep111 - ep101)/(_w2*grid[1].ticks[i1+1])) / 4
	end
	else
	( (ep010 - ep000) + (ep110 - ep100) + (ep011 - ep001) + (ep111 - ep101) ) / (4 * w2)
	end
	efv3 = ( (ep001 - ep000) + (ep101 - ep100) + (ep011 - ep010) + (ep111 - ep110) ) / (4 * w3)
	W += sum((efv1, efv2, efv3).^2) * dV * _ϵ
	end
	end
	end
	E = W * ϵ0 / 2 * u"J"
	phi_2D && (E *= 2)
	return consider_multiplicity ? E * multiplicity(grid) : E
	end

[lmh91]

Alternatively, we could drop support (prohibit) 2D-objects.

Which I am actually in favor of as our geometry is based on extended objects.

Though, for certain (rotated) objects there still could be the case where the neighbors in one dimension
of a contact grid point are both non-contact grid points, which, again, results in the 2D issue.
E.g., at a corner of a cubic geometry of a rotated contact.

[schustermartin]

Alternatively, we could drop support (prohibit) 2D-objects.

Which I am actually in favor of as our geometry is based on extended objects.

Though, for certain (rotated) objects there still could be the case where the neighbors in one dimension of a contact grid point are both non-contact grid points, which, again, results in the 2D issue. E.g., at a corner of a cubic geometry of a rotated contact.

Yes, exactly, also as the contacts are painted onto a grid, situations may occur where only grid point in a dimension is assigned as contact. Just prohibiting 2D volumes does not necessarily solve this.

[schustermartin]

I think instead of assigning electric field vectors to all grid points of the electric potential, we should assign electric field vectors to each volume between the grid points of the electric potential. Thus, the grid of the electric field would be smaller by one in each dimension as its ticks would be the midpoints of the ticks of the grid of the electric potential.

This would also require changes in the interpolation of the electric field vector during the drift.

I would not get rid of interpolation per se.
We could still interpolate, just not in voxels touching contact points. There we assign the field vector of the voxel to all points belonging to that voxel. Everywhere else interpolation can stay, and makes things smoother, especially in areas where the grid is coarse.

[lmh91]

I just looked up the Interpolations.jl doc again.
They also offer a gradient method which might be the best way to go.

using SolidStateDetectors, Unitful, Interpolations
sim = Simulation(SSD_examples[:Coax])
calculate_electric_potential!(sim)
itp = interpolate(broadcast(i -> sim.electric_potential.grid[i].ticks, (1,2,3)), sim.electric_potential.data, Gridded(Linear()));
pt = CylindricalPoint{Float32}(0.02, 0.01, 0.03)
Interpolations.gradient(itp, pt...)

Maybe we could also switch to other kind of interpolations than Linear.
At least, like @schustermartin suggested for drift steps which are not close to an contact.

We interpolate anyhow already in the drift. But we could switch to interpolating the electric field vector directly from the electric potential by this.