Strange IV curves when using PDD solver for InGaP solar cell

[UdoRoemer]

I observed an issue when using the PDD solver for InGaP solar cells. I some times get a step in the IV curve. The issue occurs and disappears (apparently) randomly when changing layer thicknesses or doping densities. In the example below, changing e.g. the base thickness to 500 nm will result in a smaller step, for some thicknesses there will be no step.
I don't get an error in the terminal and using the newest Solcore code also does not help.

So far I could not reproduce the error for other materials. I guess it is related to the known issue with high bandgap materials that is mentioned in one of the example files (there relating to the window layer, with the suggestion to remove the window from the junction).

import os
import numpy as np
import matplotlib.pyplot as plt
from solcore.light_source import LightSource
from solcore import siUnits, material, si
from solcore.solar_cell import SolarCell
from solcore.structure import Junction, Layer
from solcore.solar_cell_solver import solar_cell_solver

def this_dir_file(f):
    return os.path.join(os.path.split(__file__)[0], f)

T=298

MgF2 = material('MgF2')()
ZnS = material('ZnScub')()
InGaP = material("GaInP")

dMgF2=87.40e-9
dZnS=45.47e-9
dwindow=25e-9
demitter=100e-9
dbase=600e-9 
dbsf=100e-9 
nwindow=1e19 
nemitter=1e20 
nbase=1e17 
nbsf=1e19

V = np.linspace(-2, 2, 200)
wl = np.linspace(300, 700, 100) * 1e-9
light_source = LightSource(source_type='standard', version='AM1.5g', x=wl,
                           output_units='photon_flux_per_m', concentration=1)

user_options={'T_ambient': T, 'voltages': V, 'light_iv': True,
              'wavelength': wl, 'optics_method': 'TMM', 'light_source': light_source,
              'mpp': True,'no_back_reflection': False}



solar_cell = SolarCell(
    [
     Layer(width=dMgF2, material=MgF2, role="ARC1"),
     Layer(width=dZnS, material=ZnS, role="ARC2"),                  
     Junction([Layer(width=demitter, material=InGaP(In=0.49, Nd=siUnits(nemitter, "cm-3"), hole_diffusion_length=si("200nm"), hole_mobility = 5e-2, electron_mobility = 3.4e-3, relative_permittivity = 9), role='emitter'),
               Layer(width=dbase, material=InGaP(In=0.49, Na=siUnits(nbase, "cm-3"), electron_diffusion_length=si("200nm"), hole_mobility = 5e-2, electron_mobility = 3.4e-3, relative_permittivity = 9), role='base'),
               Layer(width=dbsf, material=InGaP(In=0.49, Na=siUnits(nbsf, "cm-3"), electron_diffusion_length=si("200nm"), hole_mobility = 5e-2, electron_mobility = 3.4e-3, relative_permittivity = 9), role='bsf')
               ], sn=100, sp=100, kind='PDD'),
    ], )

solar_cell_solver(solar_cell, 'iv', user_options)


plt.plot(V, (solar_cell.iv['IV'][1]), 'k', linewidth=4, label='Total')
print('Isc=')
print(solar_cell.iv.Isc)
print('Voc=')
print(solar_cell.iv.Voc)
print('FF=')
print(solar_cell.iv.FF)
print('Eta=')
print(solar_cell.iv.Eta)
plt.legend()
plt.ylim(-200, 200)
plt.xlim(-1.5, 1.5)
plt.ylabel('Current (A/m$^2$)')
plt.xlabel('Voltage (V)')

plt.show()


plt.semilogy(solar_cell[2].voltage, abs(solar_cell[2].current), 'k--', linewidth=4, label='Total')
plt.semilogy(solar_cell[2].voltage, abs(solar_cell[2].recombination_currents['Jrad']), 'r', label='Jrad')
plt.semilogy(solar_cell[2].voltage, abs(solar_cell[2].recombination_currents['Jsrh']), 'b', label='Jsrh')
plt.semilogy(solar_cell[2].voltage, abs(solar_cell[2].recombination_currents['Jsur']), 'g', label='Jsur')

plt.legend()
plt.xlim(-2, 2)
plt.ylim(1e-20, 1e10)
plt.xlabel('Bias (V)')
plt.ylabel('Current (A/m$^2}$)')

plt.show()

[dalonsoa]

I've run the script and got the following when solving the band structure at short circuit:

 Starting LIGHTSC...
      step  Jtot (A/m^2)          Res          Res-h    Res Poisson     Res-e        Info
         1   -0.1664E-18    0.3010E-19    0.2129E-19    0.1576E-27    0.2128E-19         2
         2   -0.1665E-17    0.3012E-18    0.2130E-18    0.1576E-27    0.2130E-18         2
         3   -0.1665E-16    0.3012E-17    0.2130E-17    0.1575E-27    0.2130E-17         2
         4   -0.1665E-15    0.3012E-16    0.2130E-16    0.1576E-27    0.2130E-16         2
         5   -0.1665E-14    0.3012E-15    0.2130E-15    0.1576E-27    0.2130E-15         3
         6   -0.1665E-13    0.3012E-14    0.2130E-14    0.1576E-27    0.2130E-14         3
         7   -0.1665E-12    0.3012E-13    0.2130E-13    0.1577E-27    0.2130E-13         3
         8   -0.1665E-11    0.3012E-12    0.2130E-12    0.1578E-27    0.2130E-12         3
         9   -0.1665E-10    0.3012E-11    0.2130E-11    0.1587E-27    0.2130E-11         3
        10   -0.1665E-09    0.3012E-10    0.2130E-10    0.1604E-27    0.2130E-10         3
        11   -0.1665E-08    0.3012E-09    0.2130E-09    0.1623E-27    0.2130E-09         3
        12   -0.1665E-07    0.3012E-08    0.2130E-08    0.1669E-27    0.2130E-08         3
        13   -0.1665E-06    0.3012E-07    0.2130E-07    0.1742E-27    0.2130E-07         3
        14   -0.1665E-05    0.3012E-06    0.2130E-06    0.1834E-27    0.2130E-06         3
        15   -0.1665E-04    0.3012E-05    0.2130E-05    0.1964E-27    0.2130E-05         3
        16   -0.1665E-03    0.3012E-04    0.2130E-04    0.1831E-27    0.2130E-04         3
        17   -0.1665E-02    0.3012E-03    0.2130E-03    0.1746E-27    0.2130E-03         3
        18   -0.1665E-01    0.3012E-02    0.2130E-02    0.1673E-27    0.2130E-02         3
        19   -0.1665        0.3012E-01    0.2130E-01    0.1623E-27    0.2130E-01         3
        20    -1.665        0.3012        0.2130        0.1605E-27    0.2130             3
        21    -16.65         3.012         2.130        0.1595E-27     2.130             3
        22    -166.5         30.12         21.30        0.1586E-27     21.30             3
 LIGHTSC Output Code: Reached Relative Tolerance.
     Res:    30.1204453
     J:   -166.515976      A/m2
     V:   -1.13121047E-32  V
 Elapsed time =   0.599416971     s

...done!

The 2nd to last column shows the residual for the transport equation for electrons, the 4th to last, the same for holes, and the last column the output code. It switches from 2 to 3 in step 4 and the residuals starts to grow. In the end there is a message saying Reached Relative Tolerance.

While reaching the relative tolerance not always means that things went wrong (it might be absolutely fine to reach relative tolerance for some problems), it indicates that things are not evolving. And a high residual is always bad news.

The summary in this case is: The problem is not converging. Anything done after this, including the calculation of the IV curve, is meaningless because the calculation of the short circuit current failed. Convergence is extremely sensitive to really small variations in the inputs. It may fail with a thickness of 500 nm and succeed if it is 501 nm. Some algorithms are more robust than others, and some might even heal themselves, going back a few steps and changing slightly the conditions. That's not the case of Solcore, unfortunately.

This is not really a bug, and there is no really a correct answer, except trying to play a bit with the inputs until you get something still compatible with the physics you are trying to model.

[dalonsoa]

In any case, it is absolutely clear that there should be:

1. A clear explanation of this in the documentation, and
2. A more evident warning of this situation, as it is done for the case of reaching the maximum number of iterations (Stop PDD solver if not convergence. Close #51 . #111 )

[UdoRoemer]

Thanks for looking into this Diego!
I am surprised that you can already see that something is off for calculations at short circuit, despite Jsc seeming to have reasonable values. So, what I get from your reply is that it will be hard to run optimisation algorithms for InGaP cells, but with some patience it might be possible to get converging simulations for some parameters.
My big question (not too much related to the original problem) is: How did you get the outputs showing the steps, Jtot and residuals?
When I solve for Isc, and look at the console, all I get is:

solar_cell_solver(solar_cell, 'short_circuit', user_options)
Solving optics of the solar cell...
Calculating RAT...
Calculating absorption profile...
Solving equilibrium...
Processing structure...
...done!

...done!

...done!

[dalonsoa]

Ok, so for your second question, I've just found out the answer: Are you executing the script as a jupyter notebook?

If so, the output of the numerical calculation (which is actually a library in fortran) will not be shown as an output cell of the notebook, but in the terminal. What the terminal is depends on how you invoked Jupyter. This is a problem because chances are that people won't look at the terminal while running a jupyter notebook, so I'll need to consider it when providing a more informative output of the PDD solver. Unfortunately, it is not straight forward to implement.

If you are not using Jupyter but you are using other IDE (eg. Spyder) it might be doing the very same thing and "swallowing" all the information from the numerical solver. None of the output you show above is actually provided by the solver, but by the Python wrapper. Try to run the script directly in a terminal with

python3 your_script.py

and let me know what you get.

About your first comment, yes, that's an issue. It is specially bad since the PDD solver (the fortran code underneath) is not designed to work in parelell, so any optmization algortihm that does it, will fail. It is a use case we didn't consider when we implemented the model back in 2014 and that is just been shown as a limitation now that more people is using Solcore and for other applications.

[UdoRoemer]

You are right, I was using Spyder. Running the script in Anaconda Prompt does the trick.
Thanks for looking into this!

[dalonsoa]

I'm closing this issue, althogh it has pointed out another issue to work on (#120)