diff --git a/.github/workflows/build_deploy_wheels.yml b/.github/workflows/build_deploy_wheels.yml index 6fb16d78..5655eee5 100644 --- a/.github/workflows/build_deploy_wheels.yml +++ b/.github/workflows/build_deploy_wheels.yml @@ -5,7 +5,7 @@ on: branches: - main - develop - - update_build_312 + tags: - "**" diff --git a/.github/workflows/test_unit_and_examples.yml b/.github/workflows/test_unit_and_examples.yml index 965d23db..0dbba4bb 100644 --- a/.github/workflows/test_unit_and_examples.yml +++ b/.github/workflows/test_unit_and_examples.yml @@ -18,7 +18,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: ["3.9", "3.10", "3.11", "3.12"] + python-version: ["3.7", "3.8", "3.9", "3.10", "3.11", "3.12"] steps: - name: Checkout @@ -53,13 +53,14 @@ jobs: - name: Install Python dependecies run: | - pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1 - python3 -m devpy install-dependencies -test-dep + pip install numpy --config-settings=setup-args="-Dallow-noblas=true" + pip install pytest meson-python ninja cython spin + python3 -m spin install-dependencies -test-dep - name: Install S4 if: matrix.os != 'windows-latest' run: | - pip install wheel + pip install wheel setuptools git clone https://github.com/phoebe-p/S4 cd S4 make S4_pyext @@ -68,14 +69,14 @@ jobs: - name: Build solcore run: | - python -m devpy build -- -Dwith_pdd=true -Dinstall_test=true + python -m spin build -- -Dwith_pdd=true -Dinstall_test=true - name: Unit and functional tests (MacOS and Linux) if: matrix.os != 'windows-latest' env: SOLCORE_SPICE: ngspice run: | - python -m devpy test -- -r a -v --cov=solcore/ --ignore=solcore/tests/test_examples.py -n "auto" + python -m spin test -- -r a -v --cov=solcore/ --ignore=solcore/tests/test_examples.py -n "auto" - name: Unit and functional tests (Windows) if: matrix.os == 'windows-latest' @@ -88,7 +89,7 @@ jobs: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} run: | python -m pip install codecov - python -m devpy codecov + python -m spin codecov test_examples: @@ -133,13 +134,13 @@ jobs: - name: Install Python dependecies run: | - pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1 - python3 -m devpy install-dependencies -test-dep + pip install pytest meson-python ninja cython numpy spin + python3 -m spin install-dependencies -test-dep - name: Install S4 if: matrix.os != 'windows-latest' run: | - pip install wheel + pip install wheel setuptools git clone https://github.com/phoebe-p/S4 cd S4 make S4_pyext @@ -148,16 +149,16 @@ jobs: - name: Build solcore run: | - python -m devpy build -- -Dwith_pdd=true -Dinstall_test=true + python -m spin build -- -Dwith_pdd=true -Dinstall_test=true - name: Unit and functional tests (MacOS and Linux) if: matrix.os != 'windows-latest' env: SOLCORE_SPICE: ngspice run: | - python -m devpy test -- -r a -v solcore/tests/test_examples.py -n "auto" + python -m spin test -- -r a -v solcore/tests/test_examples.py -n "auto" # - name: Unit and functional tests (Windows) # if: matrix.os == 'windows-latest' # run: | -# python -m devpy test -- -r a -v solcore/tests/test_examples.py +# python -m spin test -- -r a -v solcore/tests/test_examples.py diff --git a/.devpy/cmds.py b/.spin/cmds.py similarity index 96% rename from .devpy/cmds.py rename to .spin/cmds.py index 38f8153b..e63bfa1e 100644 --- a/.devpy/cmds.py +++ b/.spin/cmds.py @@ -1,6 +1,6 @@ import click -from devpy import util -from devpy.cmds.meson import _get_site_packages +from spin import util +from spin.cmds.meson import _get_site_packages @click.command() diff --git a/docs/source/Examples/DA_iv.png b/docs/source/Examples/DA_iv.png index 445d3a33..a8596d3d 100755 Binary files a/docs/source/Examples/DA_iv.png and b/docs/source/Examples/DA_iv.png differ diff --git a/docs/source/Examples/DA_qe.png b/docs/source/Examples/DA_qe.png index a3fd064f..6bc2e505 100755 Binary files a/docs/source/Examples/DA_qe.png and b/docs/source/Examples/DA_qe.png differ diff --git a/docs/source/Examples/RAT_of_ARC.png b/docs/source/Examples/RAT_of_ARC.png index c6a05278..ec6e1dfc 100755 Binary files a/docs/source/Examples/RAT_of_ARC.png and b/docs/source/Examples/RAT_of_ARC.png differ diff --git a/docs/source/Examples/example_3J_with_DA_solver.rst b/docs/source/Examples/example_3J_with_DA_solver.rst index 9661dd84..2e9fa847 100755 --- a/docs/source/Examples/example_3J_with_DA_solver.rst +++ b/docs/source/Examples/example_3J_with_DA_solver.rst @@ -6,128 +6,156 @@ Example of a 3J solar cell calculated with the DA solver .. image:: DA_qe.png :width: 40% -- Required extra files, available in `Solcore's Github repository (Examples folder) `_: - - - MgF-ZnS_AR.csv - - in01gaas.csv - - Ge-Palik.csv - .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from solcore import siUnits, material, si - from solcore.interpolate import interp1d from solcore.solar_cell import SolarCell from solcore.structure import Junction, Layer from solcore.solar_cell_solver import solar_cell_solver + from solcore.light_source import LightSource - all_materials = [] - + # Define materials for the anti-reflection coating: + MgF2 = material("MgF2")() + ZnS = material("ZnScub")() - def this_dir_file(f): - from pathlib import Path - return str(Path(__file__).parent / "data" / f) + ARC_layers = [Layer(si("100nm"), material=MgF2), + Layer(si("50nm"), material=ZnS)] + # TOP CELL - InGaP + # Now we build the top cell, which requires the n and p sides of GaInP and a window + # layer. We also add some extra parameters needed for the calculation using the + # depletion approximation: the minority carriers diffusion lengths and the doping. - # We need to build the solar cell layer by layer. - # We start from the AR coating. In this case, we load it from an an external file - refl_nm = np.loadtxt(this_dir_file("MgF-ZnS_AR.csv"), unpack=True, delimiter=",") - ref = interp1d(x=siUnits(refl_nm[0], "nm"), y=refl_nm[1], bounds_error=False, fill_value=0) - - # TOP CELL - GaInP - # Now we build the top cell, which requires the n and p sides of GaInP and a window layer. - # We also load the absorption coefficient from an external file. We also add some extra parameters needed for the - # calculation such as the minority carriers diffusion lengths AlInP = material("AlInP") InGaP = material("GaInP") window_material = AlInP(Al=0.52) - top_cell_n_material = InGaP(In=0.49, Nd=siUnits(2e18, "cm-3"), hole_diffusion_length=si("200nm")) - top_cell_p_material = InGaP(In=0.49, Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("1um")) - - all_materials.append(window_material) - all_materials.append(top_cell_n_material) - all_materials.append(top_cell_p_material) - # MID CELL - InGaAs - # We add manually the absorption coefficient of InGaAs since the one contained in the database doesn't cover - # enough range, keeping in mind that the data has to be provided as a function that takes wavelengths (m) as input and - # returns absorption (1/m) - InGaAs = material("InGaAs") - InGaAs_alpha = np.loadtxt(this_dir_file("in01gaas.csv"), unpack=True, delimiter=",") - InGaAs.alpha = interp1d(x=1240e-9 / InGaAs_alpha[0][::-1], y=InGaAs_alpha[1][::-1], bounds_error=False, fill_value=0) + top_cell_n_material = InGaP(In=0.49, Nd=siUnits(2e18, "cm-3"), + hole_diffusion_length=si("200nm")) + top_cell_p_material = InGaP(In=0.49, Na=siUnits(1e17, "cm-3"), + electron_diffusion_length=si("2um")) + # MID CELL - GaAs + GaAs = material("GaAs") - mid_cell_n_material = InGaAs(In=0.01, Nd=siUnits(3e18, "cm-3"), hole_diffusion_length=si("500nm")) - mid_cell_p_material = InGaAs(In=0.01, Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("5um")) - - all_materials.append(mid_cell_n_material) - all_materials.append(mid_cell_p_material) + mid_cell_n_material = GaAs(In=0.01, Nd=siUnits(3e18, "cm-3"), + hole_diffusion_length=si("500nm")) + mid_cell_p_material = GaAs(In=0.01, Na=siUnits(1e17, "cm-3"), + electron_diffusion_length=si("5um")) # BOTTOM CELL - Ge - # We add manually the absorption coefficient of Ge since the one contained in the database doesn't cover - # enough range. Ge = material("Ge") - Ge_alpha = np.loadtxt(this_dir_file("Ge-Palik.csv"), unpack=True, delimiter=",") - Ge.alpha = interp1d(x=1240e-9 / Ge_alpha[0][::-1], y=Ge_alpha[1][::-1], bounds_error=False, fill_value=0) - - bot_cell_n_material = Ge(Nd=siUnits(2e18, "cm-3"), hole_diffusion_length=si("800nm")) - bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("50um")) - - all_materials.append(bot_cell_n_material) - all_materials.append(bot_cell_p_material) - - # We add some other properties to the materials, assumed the same in all cases, for simplicity. - # If different, we should have added them above in the definition of the materials. - for mat in all_materials: - mat.hole_mobility = 5e-2 - mat.electron_mobility = 3.4e-3 - mat.hole_mobility = 3.4e-3 - mat.electron_mobility = 5e-2 - mat.relative_permittivity = 9 - - # And, finally, we put everything together, adding also the surface recombination velocities. We also add some shading - # due to the metallisation of the cell = 8%, and indicate it has an area of 0.7x0.7 mm2 (converted to m2) + + bot_cell_n_material = Ge(Nd=siUnits(2e18, "cm-3"), + hole_diffusion_length=si("800nm")) + bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3"), + electron_diffusion_length=si("50um")) + + # Now that the layers are configured, we can now assemble the triple junction solar + # cell. Note that we also specify a metal shading of 2% and a cell area of $1cm^{2}$. + # SolCore calculates the EQE for all three junctions and light-IV showing the relative + # contribution of each sub-cell. We set "kind = 'DA'" to use the depletion + # approximation. We can also set the surface recombination velocities, where sn + # refers to the surface recombination velocity at the n-type surface, and sp refers + # to the SRV on the p-type side. + solar_cell = SolarCell( - [ + ARC_layers + + [ Junction([Layer(si("25nm"), material=window_material, role='window'), Layer(si("100nm"), material=top_cell_n_material, role='emitter'), - Layer(si("600nm"), material=top_cell_p_material, role='base'), - ], sn=1, sp=1, kind='DA'), + Layer(si("400nm"), material=top_cell_p_material, role='base'), + ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA'), Junction([Layer(si("200nm"), material=mid_cell_n_material, role='emitter'), Layer(si("3000nm"), material=mid_cell_p_material, role='base'), - ], sn=1, sp=1, kind='DA'), + ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA'), Junction([Layer(si("400nm"), material=bot_cell_n_material, role='emitter'), Layer(si("100um"), material=bot_cell_p_material, role='base'), - ], sn=1, sp=1, kind='DA'), - ], reflectivity=ref, shading=0.08, cell_area=0.7 * 0.7 / 1e4) - - wl = np.linspace(300, 1800, 700) * 1e-9 - solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl}) + ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA') + ], + shading=0.02, cell_area=1 * 1 / 1e4) + + # Choose wavelength range (in m): + wl = np.linspace(280, 1850, 700) * 1e-9 + + # Calculate the EQE for the solar cell: + solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl, + 'da_mode': 'green', + 'optics_method': 'TMM' + }) + # we pass options to use the TMM optical method to calculate realistic R, A and T + # values with interference in the ARC (and semiconductor) layers. We can also choose + # which solver mode to use for the depletion approximation. The default is 'green', + # which uses the (faster) Green's function method. The other method is 'bvp'. + + # Plot the EQE and absorption of the individual junctions. Note that we can access + # the properties of the first junction (ignoring any other layers, such as the ARC, + # which are not part of a junction) using solar_cell(0), and the second junction using + # solar_cell(1), etc. plt.figure(1) - plt.plot(wl * 1e9, solar_cell[0].eqe(wl) * 100, 'b', label='GaInP') - plt.plot(wl * 1e9, solar_cell[1].eqe(wl) * 100, 'g', label='InGaAs') - plt.plot(wl * 1e9, solar_cell[2].eqe(wl) * 100, 'r', label='Ge') - + plt.plot(wl * 1e9, solar_cell(0).eqe(wl) * 100, 'b', label='GaInP QE') + plt.plot(wl * 1e9, solar_cell(1).eqe(wl) * 100, 'g', label='GaAs QE') + plt.plot(wl * 1e9, solar_cell(2).eqe(wl) * 100, 'r', label='Ge QE') + plt.fill_between(wl * 1e9, solar_cell(0).layer_absorption * 100, 0, alpha=0.3, + label='GaInP Abs.', color='b') + plt.fill_between(wl * 1e9, solar_cell(1).layer_absorption * 100, 0, alpha=0.3, + label='GaAs Abs.', color='g') + plt.fill_between(wl * 1e9, solar_cell(2).layer_absorption * 100, 0, alpha=0.3, + label='Ge Abs.', color='r') + + plt.plot(wl*1e9, 100*(1-solar_cell.reflected), '--k', label="100 - Reflectivity") plt.legend() plt.ylim(0, 100) plt.ylabel('EQE (%)') plt.xlabel('Wavelength (nm)') + plt.show() + + # Set up the AM0 (space) solar spectrum for the light I-V calculation: + am0 = LightSource(source_type='standard',version='AM0',x=wl, + output_units='photon_flux_per_m') + + + # Set up the voltage range for the overall cell (at which the total I-V will be + # calculated) as well as the internal voltages which are used to calculate the results + # for the individual junctions. The range of the internal_voltages should generally + # be wider than that for the voltages. - V = np.linspace(0, 3, 300) - solar_cell_solver(solar_cell, 'iv', user_options={'voltages': V, 'light_iv': True, 'wavelength': wl}) + # this is an n-p cell, so we need to scan negative voltages + + V = np.linspace(-3, 0, 300) + internal_voltages = np.linspace(-4, 2, 400) + + # Calculate the current-voltage relationship under illumination: + + solar_cell_solver(solar_cell, 'iv', user_options={'light_source': am0, + 'voltages': V, + 'internal_voltages': internal_voltages, + 'light_iv': True, + 'wavelength': wl, + 'optics_method': 'TMM', + 'mpp': True, + }) + + # We pass the same options as for solving the EQE, but also set 'light_iv' and 'mpp' to + # True to indicate we want the IV curve under illumination and to find the maximum + # power point (MPP). We also pass the AM0 light source and voltages created above. plt.figure(2) - plt.plot(V, solar_cell.iv['IV'][1], 'k', linewidth=3, label='Total') - plt.plot(V, -solar_cell[0].iv(V), 'b', label='GaInP') - plt.plot(V, -solar_cell[1].iv(V), 'g', label='InGaAs') - plt.plot(V, -solar_cell[2].iv(V), 'r', label='Ge') + plt.plot(abs(V), -solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell') + plt.plot(abs(V), solar_cell(0).iv(V)/10, 'b', label='InGaP sub-cell') + plt.plot(abs(V), solar_cell(1).iv(V)/10, 'g', label='GaAs sub-cell') + plt.plot(abs(V), solar_cell(2).iv(V)/10, 'r', label='Ge sub-cell') + plt.text(0.5,30,f'Jsc= {abs(solar_cell.iv.Isc/10):.2f} mA.cm' + r'$^{-2}$') + plt.text(0.5,28,f'Voc= {abs(solar_cell.iv.Voc):.2f} V') + plt.text(0.5,26,f'FF= {solar_cell.iv.FF*100:.2f} %') + plt.text(0.5,24,f'Eta= {solar_cell.iv.Eta*100:.2f} %') plt.legend() - plt.ylim(0, 200) + plt.ylim(0, 33) plt.xlim(0, 3) - plt.ylabel('Current (A/m$^2$)') + plt.ylabel('Current (mA/cm$^2$)') plt.xlabel('Voltage (V)') - plt.show() diff --git a/docs/source/Examples/example_RAT_of_ARC.rst b/docs/source/Examples/example_RAT_of_ARC.rst index ea2eb9ff..1d88cd8c 100755 --- a/docs/source/Examples/example_RAT_of_ARC.rst +++ b/docs/source/Examples/example_RAT_of_ARC.rst @@ -2,83 +2,69 @@ Using the TMM solver to calculate the reflection of a multilayered ARC ====================================================================== .. image:: RAT_of_ARC.png - :width: 50% - -- Required extra files, available in `Solcore's Github repository (Examples folder) `_: - - - AlInP_nk.csv - - GaInP_nk.csv - - MgF_nk.csv - - SiCx_nk.csv - - ZnS_nk.csv + :width: 100% .. code-block:: Python - from scipy.interpolate import interp1d import matplotlib.pyplot as plt import numpy as np + from solcore.solar_cell import Layer + from solcore.absorption_calculator import calculate_rat, OptiStack + from solcore import material, si - from solcore.structure import Structure - from solcore.absorption_calculator import calculate_rat - - E_eV = np.linspace(0.65, 4, 1000) - - - def nk_convert(fname, energy): - """ Designed to handle nk data files""" - - # Import data... - E_n, n, E_k, k = np.loadtxt(fname, delimiter=",", unpack=True) - print("File :: " + fname + " :: Imported Successfully!") + wl = np.linspace(300, 1900, 1000) - # Interpolate data... - n_interp = interp1d(E_n, n, bounds_error=False, fill_value=(n[0], n[-1]))(energy) - k_interp = interp1d(E_k, k, bounds_error=False, fill_value=(k[0], k[-1]))(energy) + MgF2 = material("MgF2")() + HfO2 = material("HfO2")() + ZnS = material("ZnScub")() + AlInP = material("AlInP")(Al=0.52) + GaInP = material("GaInP")(In=0.49) - return (energy, n_interp, k_interp) - - - # Load nk data using nk_convert function... - alinp_nk = nk_convert("AlInP_nk.csv", energy=E_eV) - gainp_nk = nk_convert("GaInP_nk.csv", energy=E_eV) - mgf_nk = nk_convert("MgF_nk.csv", energy=E_eV) - sic_nk = nk_convert("SiCx_nk.csv", energy=E_eV) - zns_nk = nk_convert("ZnS_nk.csv", energy=E_eV) - - # Build the optical stack... - stack = Structure([ - [117, 1240 / E_eV, mgf_nk[1], mgf_nk[2]], - [80, 1240 / E_eV, sic_nk[1], sic_nk[2]], - [61, 1240 / E_eV, zns_nk[1], zns_nk[2]], - [25, 1240 / E_eV, alinp_nk[1], alinp_nk[2]], - [350000, 1240 / E_eV, gainp_nk[1], gainp_nk[2]] - ]) + stack = OptiStack([ + Layer(si('141nm'), material=MgF2), + Layer(si('82nm'), material=HfO2), + Layer(si('70nm'), material=ZnS), + Layer(si('25nm'), material=AlInP), + ], substrate=GaInP, no_back_reflection=False) angles = np.linspace(0, 80, 10) - RAT_angles = [] + RAT_angles = np.zeros((len(angles), 3, len(wl))) - print("Calculate RAT ::") - for theta in angles: - print("Calculating at angle :: %4.1f deg" % theta) + print("Calculate RAT:") + for i, theta in enumerate(angles): + print("Calculating at angle: %4.1f deg" % theta) # Calculate RAT data... - rat_data = calculate_rat(stack, angle=theta, wavelength=1240 / E_eV) + rat_data = calculate_rat(stack, angle=theta, wavelength=wl, + no_back_reflection=False) - RAT_angles.append((theta, rat_data["R"], rat_data["A"])) + RAT_angles[i] = [rat_data["R"], rat_data["A"], rat_data["T"]] colors = plt.cm.jet(np.linspace(1, 0, len(RAT_angles))) - fig, ax2 = plt.subplots(1, 1) + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) for i, RAT in enumerate(RAT_angles): - ax2.plot(1240 / E_eV, RAT[1] * 100, ls="-", color=colors[i], label="%4.1f$^\circ$" % RAT[0]) - ax2.plot(1240 / E_eV, RAT[2] * 100, ls="--", color=colors[i]) + if i == 0: + ax1.plot(wl, RAT[0] * 100, ls="-", color=colors[i], label="R") + ax1.plot(wl, (RAT[1] + RAT[2]) * 100, ls="--", color=colors[i], label="A+T") + + else: + ax1.plot(wl, RAT[0] * 100, ls="-", color=colors[i]) + ax1.plot(wl, (RAT[1] + RAT[2]) * 100, ls="--", color=colors[i]) + + ax2.plot(wl, RAT[1]*100, color=colors[i], label="%4.1f$^\circ$" % angles[i]) + + ax1.set_ylim([0, 100]) + ax1.set_xlim([300, 1800]) + ax1.set_xlabel("Wavelength (nm)") + ax1.set_ylabel("Reflection and transmission into structure (%)") + ax1.legend(loc=5) ax2.set_ylim([0, 100]) ax2.set_xlim([300, 1800]) ax2.set_xlabel("Wavelength (nm)") - ax2.set_ylabel("Reflection and Transmission (%)") + ax2.set_ylabel("Absorption in surface layers (%)") ax2.legend(loc=5) - ax2.text(0.05, 0.45, '(a)', transform=ax2.transAxes, fontsize=12) - plt.tight_layout(w_pad=4) - plt.show() + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/docs/source/Examples/main.rst b/docs/source/Examples/main.rst index 94ffea4b..8e8ced00 100755 --- a/docs/source/Examples/main.rst +++ b/docs/source/Examples/main.rst @@ -17,7 +17,7 @@ Light sources Optical methods --------------- -- :doc:`Using the TMM solver to calculate the reflextion of a multilayered ARC ` +- :doc:`Using the TMM solver to calculate the reflection of a multilayered ARC ` - :doc:`Looking at the effect of substrate and the no_back_reflection option in the TMM solver ` - :doc:`Comparing optical models (TMM, Beer-Lambert, and RCWA) ` diff --git a/docs/source/Solvers/DDsolver.rst b/docs/source/Solvers/DDsolver.rst index d3267778..1528fb8b 100755 --- a/docs/source/Solvers/DDsolver.rst +++ b/docs/source/Solvers/DDsolver.rst @@ -3,12 +3,20 @@ Fortran Poisson Drift-Diffusion solver (PDD) This section documents the functionality of the original, Fortran-based Poisson Drift-Diffusion Solver included in Solcore. For the newer interface to the Python-based Sesame solver, which also solves the PDD equations, see -:doc:`sesame`. +:doc:`sesame`. While the two solvers should in principle have the same functionality, the Fortran solver was developed +specifically to deal with solar cells containing quantum wells (QWs), where the abrupt and frequent change in material parameters +can cause convergence issues. The Sesame solver has not been tested for the QW use case. The Fortran solver is less +suitable for cells containing thick (> 10 microns or so) layers, as it was designed for III-V cells. The Sesame solver +was designed for silicon-based cells and thus can handle thicker layers. The Fortran solver can handle only constant +values of doping per layer, while the Sesame solver also accepts depth-dependent doping profiles defined by the user. - Example 1: :doc:`Example of a simple 2J solar cell calculated with the PDD solver <../Examples/example_PDD_solver>` - Example 2: :doc:`Tutorial: 2J solar cell with QWs in the bottom cell <../Examples/tutorial>` -The PDD package provide all tools necesary to build a solar cell structure and calculate its properties by solving simultaneously the Poisson equation and the drfit diffusion equations. Normally, these functions will not need to be accessed directly, but are called internally by :literal:`Solcore` when using the higher level methods in the :doc:`solar cell solver `. +The PDD package provide all tools necesary to build a solar cell structure and calculate its properties by solving +simultaneously the Poisson equation and the drfit diffusion equations. Normally, these functions will not need to be +accessed directly, but are called internally by :literal:`Solcore` when using the higher level methods in the +:doc:`solar cell solver `. For using the PDD package, it is enough to include the following line in your code: diff --git a/docs/source/Solvers/Figures/PDD_np.png b/docs/source/Solvers/Figures/PDD_np.png new file mode 100644 index 00000000..904a6089 Binary files /dev/null and b/docs/source/Solvers/Figures/PDD_np.png differ diff --git a/docs/source/Solvers/Figures/PDD_pn.png b/docs/source/Solvers/Figures/PDD_pn.png new file mode 100644 index 00000000..ee030112 Binary files /dev/null and b/docs/source/Solvers/Figures/PDD_pn.png differ diff --git a/docs/source/Solvers/SesameDDsolver.rst b/docs/source/Solvers/SesameDDsolver.rst index 434c259f..f39682b4 100644 --- a/docs/source/Solvers/SesameDDsolver.rst +++ b/docs/source/Solvers/SesameDDsolver.rst @@ -1,3 +1,91 @@ -Poisson - Drift-Diffusion solver (PDD) +Sesame Poisson Drift-Diffusion solver ====================================== +This solver provides an interface to `Sesame `_ to solve the +Drift-Diffusion equations. While it serves the same purpose as the legacy Fortran PDD solver, this solver has +the advantage of being written entirely in Python, thus being more transparent to most users. +While the two solvers should in principle have the same functionality, the Fortran solver was developed +specifically to deal with solar cells containing quantum wells (QWs), where the abrupt and frequent change in material parameters +can cause convergence issues. The Sesame solver has not been tested for the QW use case. The Fortran solver is less +suitable for cells containing thick (> 10 microns or so) layers, as it was designed for III-V cells. The Sesame solver +was designed for silicon-based cells and thus can handle thicker layers. + +Unlike the Fortran-based solver, this solver can handle both constant values of doping per layer, +or depth-dependent doping profiles defined by the user. + +Material constants +------------------- + +The following material parameters will be extracted from each layer in a junction and used by +Sesame: + +- ``Nc``: Conduction band effective density of states +- ``Nv``: Valence band effective density of states +- ``band_gap`` +- ``electron_affinity`` +- ``relative_permittivity`` +- ``electron_mobility`` +- ``hole_mobility`` +- ``electron_minority_lifetime`` +- ``hole_minority_lifetime`` +- ``bulk_recombination_energy``: this is an optional parameter (will be set to 0 by default). + It should be set (in joules!), in the ``material`` definition, if the user wants to use a different value +- ``radiative_recombination``: radiative recombination rate (m\ :sup:`3` s\ :sup:`-1`) +- ``electron_auger_recombination``: Auger recombination rate for electrons (m\ :sup:`6` s\ :sup:`-1`) +- ``hole_auger_recombination``: Auger recombination rate for holes (m\ :sup:`6` s\ :sup:`-1`) + +Note that this list uses the names of the parameters used by Solcore's material system, not the +names used internally by Sesame (which can be found in the Sesame documentation). While Sesame internally uses different units, values should be passed to +materials/layers in base SI units, as for the rest of Solcore (m, s, kg, etc). + +Mesh +----- + +Solcore will try to construct a reasonable mesh to solve the PDD equations, putting +more mesh points close to the front surface of the cell and in regions where the doping +is changing. However, the user can also provide a custom mesh, defined in terms of +distance *from the front surface of the junction* in m. +This is is passed as a 1-dimensional array through the ``mesh`` argument to ``Junction``). + +Doping profile +--------------- + +As for the Fortran PDD solver and the depletion approximation (DA) solver, the user can set fixed +doping levels for each layer in a junction using the ``Nd`` argument for n-type doping and the ``Na`` +argument for p-type. However, it is also possible to define depth-dependent doping profiles. These should be a function +which accepts an array of positions (in m) and returns the doping at that position in units of +m\ :sup:`-3`. Sesame will interpret positive values as n-type doping and negative values as p-type. +The doping profile can be specified for the whole junction by setting the ``doping_profile`` argument +of ``Junction``, or for individual layers by setting the ``doping_profile`` argument for ``Layer``. It is +possible to mix constant doping in some layers with depth-dependent doping in others. The position argument +of the function should always be in terms of distance from the front surface of the junction or layer it is for. + +Outputs +------- +In addition to updating the individual junctions and the solar cell as a whole with the current-voltage output +(if ``solar_cell_solver`` is called with task ``iv``) and quantum efficiency +(if ``solar_cell_solver`` is called with task ``qe``), the Sesame solver will also update each junction it solves the +IV for with an attribute called ``pdd_output``, which contains the following: + +- ``pdd_output.potential``: the potential (V) +- ``pdd_output.n``: the electron density (m\ :sup:`-3`) +- ``pdd_output.p``: the hole density (m\ :sup:`-3`) +- ``pdd_output.Ec``: the level of the conduction band (eV) +- ``pdd_output.Ev``: the level of the valence band (eV) +- ``pdd_output.Efe``: the electron quasi-Fermi level (eV) +- ``pdd_output.Efh``: the hole quasi-Fermi level (eV) +- ``pdd_output.Rrad``: the radiative recombination rate (m\ :sup:`-3` s\ :sup:`-1`) +- ``pdd_output.Raug``: the Auger recombination rate (m\ :sup:`-3` s\ :sup:`-1`) +- ``pdd_output.Rsrh``: the bulk recombination due to Shockley-Read-Hall processes (m\ :sup:`-3` s\ :sup:`-1`) + +Each of these is a 2-dimensional array, with dimensions ``(len(options.internal_voltages), len(mesh))``. + +Sub-function documentation +--------------------------- + +Note that the functions below should generally not be accessed directly, but will be called by +``solar_cell_solver`` for junctions with ``kind="sesame_PDD"``. + +.. automodule:: solcore.sesame_drift_diffusion.solve_pdd + :members: + :undoc-members: \ No newline at end of file diff --git a/docs/source/Solvers/sign_conventions.rst b/docs/source/Solvers/sign_conventions.rst index 63ab6367..538e3020 100644 --- a/docs/source/Solvers/sign_conventions.rst +++ b/docs/source/Solvers/sign_conventions.rst @@ -1,9 +1,33 @@ +.. _sign-conventions: + Current-voltage sign conventions ================================ In order to maintain internal consistency, the following sign conventions are used for current and voltage: - A positive voltage means that a higher potential is applied to the front contact of the cell. - For a p-n junction (i.e. p-type material is on top of n-type material), this means that a positive - voltage results in 'solar cell operation' under illumination, with a positive Voc. - For an n-p junction, a negative voltage must be applied and Voc will be negative. + + - For a **p-n junction** (i.e. p-type material is on top of n-type material), this means that a **positive + voltage** results in 'solar cell operation' under illumination, with a positive open-circuit voltage (V\ :sub:`OC`). + - For an **n-p junction**, a **negative voltage** must be applied for PV operation and V\ :sub:`OC` will be negative. + +- For the junction object itself, the short-circuit current (J\ :sub:`SC`) is negative if the cell operates at positive voltage, and positive + if the cell operates at negative voltage (the product of current and voltage is negative, i.e. work is being done + *by* the cell). This relates to the data which is stored in ``junction.iv``. + +- For the cell as a whole (i.e. the multi-junction IV for a multi-junction cell, or the IV with shunt resistance for + a single-junction cell), the short-circuit current was taken to be positive for positive voltages. This means that + it will be negative for negative voltages. + +- The two-diode (2D) and detailed-balance (DB) models always assume that, in this sign convention, the cell has the polarity of + a p-n junction (positive voltage should be applied) + +.. image:: Figures/PDD_pn.png + :width: 40% +.. image:: Figures/PDD_np.png + :width: 40% + +Note that the sign conventions here were always the case for the (Fortran) PDD solver, but were not used consistently +by the depletion approximation solver before version 5.10, which always assumed that V\ :sub:`OC` and J\ :sub:`SC` were positive. This +does not matter for single-junction cells, since the choice is arbitrary, but caused inconsistencies when defining +n-p cells which mix different types of junction. diff --git a/docs/source/Solvers/solving_solar_cells.rst b/docs/source/Solvers/solving_solar_cells.rst index 898a50d8..6b7c5c00 100755 --- a/docs/source/Solvers/solving_solar_cells.rst +++ b/docs/source/Solvers/solving_solar_cells.rst @@ -21,18 +21,33 @@ The *task* has to be "optics", "iv", "qe", "equilibrium" or "short_circuit", the All of the above calculations modify the original solar cell object, adding attributes or methods to its structure. For example, after calculating the IV curve of my_solar_cell_object, this will have a new attribute called "iv" that is a dictionary with the IV of the total solar cell, the IV curves of each of the junctions, information of the Voc and Isc, if relevant, etc. -More details of the specific electrical solvers included in Solcore can be found in: +Kinds of junction +----------------- + +The following types of junctions are currently implemented in Solcore, specified through the :literal:`kind` argument passed to +when creating a :literal:`Junction` object: + +1. **DB**: detailed balance (the Shockley-Queisser limit) +2. **2D**: two-diode model +3. **DA**: depletion approximation +4. **PDD**: Poisson-drift-diffusion, using the legacy Fortran solver +5. **sesame_PDD**: Poisson-drift-diffusion, using the new Python-based Sesame solver + +For the current and voltage sign conventions used by Solcore, please see :doc:`here `. + +More details of these specific electrical solvers included in Solcore can be found in: .. toctree:: :maxdepth: 0 + sign_conventions detailed_balance TwoDiode depletion + SesameDDsolver DDsolver multijunction_iv - .. _solver-options: Solver Options diff --git a/examples/GaAs_cell_drift_diffusion.py b/examples/GaAs_cell_drift_diffusion.py deleted file mode 100644 index d6fe102f..00000000 --- a/examples/GaAs_cell_drift_diffusion.py +++ /dev/null @@ -1,162 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from solcore import material, si -from solcore.structure import Junction, Layer -from solcore.solar_cell_solver import solar_cell_solver, SolarCell -from solcore.state import State -from solcore.light_source import LightSource - -# set up structure similar to https://doi.org/10.1109/16.46385 - -# define user options -options = State() -options.wavelength = np.linspace(280, 950, 100)*1e-9 -options.voltages = np.linspace(0, 1.3, 20) -options.light_iv = True -options.mpp = True -options.optics_method = 'TMM' -options.recalculate_absorption = True -options.light_source = LightSource(source_type="standard", - version="AM1.5g", x=options.wavelength, output_units="photon_flux_per_m") - - -# define materials - -MgF2 = material("MgF2")() -ZnS = material("ZnScub")() - -ARC_layers = [Layer(si("100nm"), material=MgF2), - Layer(si("50nm"), material=ZnS)] - -AlGaAs_window = material('AlGaAs')(Al=0.8, Na=si('4e18cm-3'), - relative_permittivity=11, - electron_minority_lifetime=1e-6, - hole_minority_lifetime=1e-6 - ) -GaAs_emitter = material('GaAs')(Na=si('4e18cm-3'), - electron_minority_lifetime=1e-6, - hole_minority_lifetime=1e-6) -GaAs_base = material('GaAs')(Nd=si('2e17cm-3'), - electron_minority_lifetime=1e-6, - hole_minority_lifetime=1e-6) -GaAs_bsf = material('GaAs')(Nd=si('2e18cm-3'), - electron_minority_lifetime=1e-6, - hole_minority_lifetime=1e-6 - ) -GaAs_substrate = material('GaAs')() - -junction_layers = [ - Layer(si('30nm'), AlGaAs_window, role='window'), - Layer(si('100nm'), GaAs_emitter, role='emitter'), - Layer(si('3000nm'), GaAs_base, role='base'), - Layer(si('100nm'), GaAs_bsf, role='bsf'), - ] - -# define solar cell -mesh = np.linspace(0, 3230, 4000)*1e-9 - -solar_cell = SolarCell( - ARC_layers + [Junction(junction_layers, kind='PDD', sn=1e6, sp=1e6,)], - substrate=GaAs_substrate) - -solar_cell_sesame = SolarCell( - ARC_layers + [Junction(junction_layers, kind='sesame_PDD', sn=1e6, sp=1e6, - # mesh=mesh - )], - substrate=GaAs_substrate) - -solar_cell_optics = SolarCell( - ARC_layers + junction_layers, substrate=GaAs_substrate -) - -solar_cell_solver(solar_cell_optics, 'optics', options) - -solar_cell_solver(solar_cell, 'iv', options) -solar_cell_solver(solar_cell, 'qe', options) - -solar_cell_solver(solar_cell_sesame, 'iv', options) -solar_cell_solver(solar_cell_sesame, 'qe', options) - -absorption_per_layer = np.array([layer.layer_absorption for layer in solar_cell_optics]) - -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) - -ax1.plot(solar_cell.iv['IV'][0], solar_cell.iv['IV'][1]/10, label='IV') -ax1.plot(options.voltages, solar_cell[2].iv(options.voltages)/10, label='IV internal') - -ax2.stackplot(options.wavelength*1e9, 100*absorption_per_layer[::-1], alpha=0.4) -ax2.plot(options.wavelength*1e9, 100*solar_cell[2].eqe(options.wavelength), '-k') - -ax2.legend(['GaAs BSF', 'GaAs base', 'GaAs emitter', 'AlGaAs window', 'ZnS', 'MgF$_2$']) - -ax1.set_xlabel('Voltage (V)') -ax1.set_ylabel('Current density (mA/cm$^2$)') -ax1.set_ylim(-10, 35) - -ax2.set_xlabel('Wavelength (nm)') -ax2.set_ylabel('EQE (%)') - -plt.tight_layout() -plt.show() - -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) - -ax1.plot(solar_cell_sesame.iv['IV'][0], solar_cell_sesame.iv['IV'][1]/10, label='IV') -ax1.plot(options.voltages, solar_cell_sesame[2].iv(options.voltages)/10, label='IV internal') - -ax2.stackplot(options.wavelength*1e9, 100*absorption_per_layer[::-1], alpha=0.4) -ax2.plot(options.wavelength*1e9, 100*solar_cell_sesame[2].eqe(options.wavelength), '-k') - -ax2.legend(['GaAs BSF', 'GaAs base', 'GaAs emitter', 'AlGaAs window', 'ZnS', 'MgF$_2$']) - -ax1.set_xlabel('Voltage (V)') -ax1.set_ylabel('Current density (mA/cm$^2$)') -ax1.set_ylim(-10, 35) - -ax2.set_xlabel('Wavelength (nm)') -ax2.set_ylabel('EQE / Absorption (%)') - -plt.tight_layout() -plt.show() - -# from solcore import material, si -# from solcore.structure import Junction, Layer -# from solcore.solar_cell_solver import solar_cell_solver, SolarCell -# from solcore.analytic_solar_cells import iv_depletion -# from solcore.sesame_drift_diffusion import iv_sesame -# from solcore.state import State -# from solcore.optics import solve_tmm -# from solcore.light_source import LightSource -# import numpy as np -# -# options = State() -# options.wavelength = np.linspace(300, 950, 100) * 1e-9 -# options.voltages = np.linspace(-1.2, 1.2, 20) -# options.internal_voltages = np.linspace(-1.2, 1.2, 20) -# options.T = 300 -# options.light_iv = True -# options.light_source = LightSource(source_type='standard', version='AM1.5g', x=options.wavelength, -# output_units='photon_flux_per_m') -# options.da_mode = 'green' -# options.optics_method = 'TMM' -# -# mesh = np.linspace(0, 3100, 2000)*1e-9 -# -# GaAs_n = material('GaAs')(T=300, Na=si('4e18cm-3'), hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6, -# electron_mobility=1, hole_mobility=1) -# GaAs_p = material('GaAs')(T=300, Nd=si('2e17cm-3'), hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6, -# electron_mobility=1, hole_mobility=1) -# -# pn_junction = Junction([Layer(si('100nm'), GaAs_n, role='emitter'), Layer(si('3000nm'), GaAs_p, role='base')], -# kind='sesame_PDD', -# # mesh=mesh -# ) -# -# solar_cell_solver(SolarCell([pn_junction]), 'optics', options) -# -# plt.figure() -# iv_sesame(pn_junction, options) -# plt.plot(options.voltages, pn_junction.iv(options.voltages), label='Sesame') -# plt.ylim(-250, 250) -# plt.show() - diff --git a/examples/GaAs_cell_drift_diffusion_np.py b/examples/GaAs_cell_drift_diffusion_np.py new file mode 100644 index 00000000..584ad08a --- /dev/null +++ b/examples/GaAs_cell_drift_diffusion_np.py @@ -0,0 +1,137 @@ +import numpy as np +import matplotlib.pyplot as plt +from solcore import material, si +from solcore.structure import Junction, Layer +from solcore.solar_cell_solver import solar_cell_solver, SolarCell +from solcore.state import State +from solcore.light_source import LightSource + +# set up structure similar to https://doi.org/10.1109/16.46385 (but with doping flipped) +# Compare performance of the two PDD solvers. + +# define user options +options = State() +options.wavelength = np.linspace(280, 950, 70)*1e-9 +options.voltages = np.linspace(-1.2, -0.2, 60) +options.internal_voltages = np.linspace(-1.2, 0.00, 60) +options.light_iv = True +options.mpp = True +options.optics_method = 'TMM' +options.light_source = LightSource(source_type="standard", + version="AM1.5g", x=options.wavelength, output_units="photon_flux_per_m") +options.no_back_reflection = False +options.recalculate_absorption = True + + +# define materials + +MgF2 = material("MgF2")() +ZnS = material("ZnScub")() +Ag = material("Ag")() + +ARC_layers = [Layer(si("100nm"), material=MgF2), + Layer(si("50nm"), material=ZnS)] + +AlGaAs_window = material('AlGaAs')(Al=0.8, Nd=si('4e18cm-3'), + relative_permittivity=11, + electron_minority_lifetime=1e-9, + hole_minority_lifetime=1e-9, + ) +GaAs_emitter = material('GaAs')(Nd=si('4e18cm-3'), + electron_minority_lifetime=5e-9, + hole_minority_lifetime=5e-9) +GaAs_base = material('GaAs')(Na=si('2e17cm-3'), + electron_minority_lifetime=5e-9, + hole_minority_lifetime=5e-9) +GaAs_bsf = material('GaAs')(Na=si('2e18cm-3'), + electron_minority_lifetime=5e-9, + hole_minority_lifetime=5e-9, + ) +GaAs_substrate = material('GaAs')() + +junction_layers = [ + Layer(si('30nm'), AlGaAs_window, role='window'), + Layer(si('100nm'), GaAs_emitter, role='emitter'), + Layer(si('3000nm'), GaAs_base, role='base'), + Layer(si('100nm'), GaAs_bsf, role='bsf'), + ] + +solar_cell = SolarCell( + ARC_layers + [Junction(junction_layers, kind='PDD', sn=1e6, sp=1e6, + R_shunt=1, + ), + Layer(si('100um'), GaAs_substrate) + ], + substrate=Ag, + R_series=0.0001, +) + +solar_cell_sesame = SolarCell( + ARC_layers + [Junction(junction_layers, kind='sesame_PDD', sn=1e6, sp=1e6, + R_shunt=1, + ), + Layer(si('100um'), GaAs_substrate) + ], + substrate=Ag, + R_series=0.0001, +) + +solar_cell_optics = SolarCell( + ARC_layers + junction_layers + [Layer(si('100um'), GaAs_substrate)], substrate=Ag, +) + +solar_cell_solver(solar_cell, 'iv', options) +solar_cell_solver(solar_cell, 'qe', options) + +solar_cell_solver(solar_cell_sesame, 'iv', options) +solar_cell_solver(solar_cell_sesame, 'qe', options) + +solar_cell_solver(solar_cell_optics, 'optics', options) + +absorption_per_layer = np.array([layer.layer_absorption for layer in solar_cell_optics]) + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) + +ax1.plot(solar_cell.iv['IV'][0], solar_cell.iv['IV'][1]/10, '--k', label='Overall IV') +ax1.plot(options.voltages, -solar_cell[2].iv(options.voltages)/10, '-r', label='junction IV') + +ax1.legend() +ax1.set_title("Fortran PDD solver") + +ax2.stackplot(options.wavelength*1e9, 100*absorption_per_layer[::-1], alpha=0.4) +ax2.plot(options.wavelength*1e9, 100*solar_cell[2].eqe(options.wavelength), '-k') + +ax2.legend(['substrate', 'GaAs BSF', 'GaAs base', 'GaAs emitter', 'AlGaAs window', 'ZnS', 'MgF$_2$']) + +ax1.set_xlabel('Voltage (V)') +ax1.set_ylabel('Current density (mA/cm$^2$)') +ax1.set_ylim(-35, 10) + +ax2.set_xlabel('Wavelength (nm)') +ax2.set_ylabel('EQE / Absorption (%)') + +plt.tight_layout() +plt.show() + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) + +ax1.plot(solar_cell_sesame.iv['IV'][0], solar_cell_sesame.iv['IV'][1]/10, '--k', label='Overall IV') +ax1.plot(options.voltages, -solar_cell_sesame[2].iv(options.voltages)/10, '-r', label='junction IV') + +ax1.legend() +ax1.set_title("Sesame PDD solver") + +ax2.stackplot(options.wavelength*1e9, 100*absorption_per_layer[::-1], alpha=0.4) +ax2.plot(options.wavelength*1e9, 100*solar_cell_sesame[2].eqe(options.wavelength), '-k') + +ax2.legend(['substrate', 'GaAs BSF', 'GaAs base', 'GaAs emitter', 'AlGaAs window', 'ZnS', 'MgF$_2$']) + +ax1.set_xlabel('Voltage (V)') +ax1.set_ylabel('Current density (mA/cm$^2$)') +ax1.set_ylim(-35, 10) + +ax2.set_xlabel('Wavelength (nm)') +ax2.set_ylabel('EQE / Absorption (%)') + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/examples/GaAs_cell_drift_diffusion_pn.py b/examples/GaAs_cell_drift_diffusion_pn.py new file mode 100644 index 00000000..718842ac --- /dev/null +++ b/examples/GaAs_cell_drift_diffusion_pn.py @@ -0,0 +1,135 @@ +import numpy as np +import matplotlib.pyplot as plt +from solcore import material, si +from solcore.structure import Junction, Layer +from solcore.solar_cell_solver import solar_cell_solver, SolarCell +from solcore.state import State +from solcore.light_source import LightSource + +# set up structure similar to https://doi.org/10.1109/16.46385 + +# define user options +options = State() +options.wavelength = np.linspace(280, 950, 70)*1e-9 +options.voltages = np.linspace(0, 1.1, 30) +options.internal_voltages = options.voltages +options.light_iv = True +options.mpp = True +options.optics_method = 'TMM' +options.light_source = LightSource(source_type="standard", + version="AM1.5g", x=options.wavelength, output_units="photon_flux_per_m") +options.no_back_reflection = False +options.recalculate_absorption = True + +# define materials + +MgF2 = material("MgF2")() +ZnS = material("ZnScub")() +Ag = material("Ag")() + +ARC_layers = [Layer(si("100nm"), material=MgF2), + Layer(si("50nm"), material=ZnS)] + +AlGaAs_window = material('AlGaAs')(Al=0.8, Na=si('4e18cm-3'), + relative_permittivity=11, + electron_minority_lifetime=1e-9, + hole_minority_lifetime=1e-9, + ) +GaAs_emitter = material('GaAs')(Na=si('4e18cm-3'), + electron_minority_lifetime=5e-9, + hole_minority_lifetime=5e-9) +GaAs_base = material('GaAs')(Nd=si('2e17cm-3'), + electron_minority_lifetime=5e-9, + hole_minority_lifetime=5e-9) +GaAs_bsf = material('GaAs')(Nd=si('2e18cm-3'), + electron_minority_lifetime=5e-9, + hole_minority_lifetime=5e-9, + ) +GaAs_substrate = material('GaAs')() + +junction_layers = [ + Layer(si('30nm'), AlGaAs_window, role='window'), + Layer(si('100nm'), GaAs_emitter, role='emitter'), + Layer(si('3000nm'), GaAs_base, role='base'), + Layer(si('100nm'), GaAs_bsf, role='bsf'), + ] + +solar_cell = SolarCell( + ARC_layers + [Junction(junction_layers, kind='PDD', sn=1e6, sp=1e6, + R_shunt=1, + ), + Layer(si('100um'), GaAs_substrate) + ], + substrate=Ag, + R_series=0.0001, +) + +solar_cell_sesame = SolarCell( + ARC_layers + [Junction(junction_layers, kind='sesame_PDD', sn=1e6, sp=1e6, + R_shunt=1, + ), + Layer(si('100um'), GaAs_substrate) + ], + substrate=Ag, + R_series=0.0001, +) + +solar_cell_optics = SolarCell( + ARC_layers + junction_layers + [Layer(si('100um'), GaAs_substrate)], substrate=Ag, +) + +solar_cell_solver(solar_cell, 'iv', options) +solar_cell_solver(solar_cell, 'qe', options) + +solar_cell_solver(solar_cell_sesame, 'iv', options) +solar_cell_solver(solar_cell_sesame, 'qe', options) + +solar_cell_solver(solar_cell_optics, 'optics', options) + +absorption_per_layer = np.array([layer.layer_absorption for layer in solar_cell_optics]) + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) + +ax1.plot(solar_cell.iv['IV'][0], solar_cell.iv['IV'][1]/10, '--k', label='Overall IV') +ax1.plot(options.voltages, -solar_cell[2].iv(options.voltages)/10, '-r', label='junction IV') + +ax1.legend() +ax1.set_title("Fortran PDD solver") + +ax2.stackplot(options.wavelength*1e9, 100*absorption_per_layer[::-1], alpha=0.4) +ax2.plot(options.wavelength*1e9, 100*solar_cell[2].eqe(options.wavelength), '-k') + +ax2.legend(['substrate', 'GaAs BSF', 'GaAs base', 'GaAs emitter', 'AlGaAs window', 'ZnS', 'MgF$_2$']) + +ax1.set_xlabel('Voltage (V)') +ax1.set_ylabel('Current density (mA/cm$^2$)') +ax1.set_ylim(-10, 35) + +ax2.set_xlabel('Wavelength (nm)') +ax2.set_ylabel('EQE / Absorption (%)') + +plt.tight_layout() +plt.show() + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) + +ax1.plot(solar_cell_sesame.iv['IV'][0], solar_cell_sesame.iv['IV'][1]/10, '--k', label='Overall IV') +ax1.plot(options.voltages, -solar_cell_sesame[2].iv(options.voltages)/10, '-r', label='junction IV') + +ax1.legend() +ax1.set_title("Sesame PDD solver") + +ax2.stackplot(options.wavelength*1e9, 100*absorption_per_layer[::-1], alpha=0.4) +ax2.plot(options.wavelength*1e9, 100*solar_cell_sesame[2].eqe(options.wavelength), '-k') + +ax2.legend(['substrate', 'GaAs BSF', 'GaAs base', 'GaAs emitter', 'AlGaAs window', 'ZnS', 'MgF$_2$']) + +ax1.set_xlabel('Voltage (V)') +ax1.set_ylabel('Current density (mA/cm$^2$)') +ax1.set_ylim(-10, 35) + +ax2.set_xlabel('Wavelength (nm)') +ax2.set_ylabel('EQE / Absorption (%)') + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/examples/MJ_solar_cell_PDD_solver.py b/examples/MJ_solar_cell_PDD_solver.py index be5f2875..d1b0dab4 100755 --- a/examples/MJ_solar_cell_PDD_solver.py +++ b/examples/MJ_solar_cell_PDD_solver.py @@ -16,41 +16,45 @@ # TOP CELL - InGaP # Now we build the top cell, which requires the n and p sides of GaInP and a window -# layer. We also add some extra parameters needed for the calculation using the -# depletion approximation: the minority carriers diffusion lengths and the doping. +# layer. We can also add extra parameters (in this case, the carrier lifetimes and +# the doping levels). AlInP = material("AlInP") InGaP = material("GaInP") -window_material = AlInP(Al=0.52, Nd=siUnits(2e18, "cm-3"), relative_permittivity=9, - electron_minority_lifetime=1e-9, hole_minority_lifetime=1e-9) +window_material = AlInP(Al=0.52, Nd=siUnits(3e18, "cm-3"), relative_permittivity=9, + electron_minority_lifetime=1e-10, hole_minority_lifetime=1e-10) top_cell_n_material = InGaP(In=0.49, Nd=siUnits(2e18, "cm-3"), electron_minority_lifetime=5e-9, hole_minority_lifetime=5e-9) top_cell_p_material = InGaP(In=0.49, Na=siUnits(1e17, "cm-3"), electron_minority_lifetime=5e-9, hole_minority_lifetime=5e-9) -# For convenience we will set the carrier transport properties to default values for -# all the materials in one place. To do that we add the materials to a list called -# all_materials: - # MID CELL - GaAs GaAs = material("GaAs") -mid_cell_n_material = GaAs(In=0.01, Nd=siUnits(3e18, "cm-3")) -mid_cell_p_material = GaAs(In=0.01, Na=siUnits(1e17, "cm-3")) +mid_cell_n_material = GaAs(In=0.01, Nd=siUnits(3e18, "cm-3"), + electron_minority_lifetime=5e-9, hole_minority_lifetime=5e-9) +mid_cell_p_material = GaAs(In=0.01, Na=siUnits(1e17, "cm-3"), + electron_minority_lifetime=5e-9, hole_minority_lifetime=5e-9) # BOTTOM CELL - Ge Ge = material("Ge") -bot_cell_n_material = Ge(Nd=siUnits(2e18, "cm-3")) -bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3")) - +bot_cell_n_material = Ge(Nd=siUnits(2e18, "cm-3"), + electron_minority_lifetime=5e-5, hole_minority_lifetime=5e-5) +bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3"), + electron_minority_lifetime=5e-4, hole_minority_lifetime=5e-4) # Now that the layers are configured, we can now assemble the triple junction solar # cell. Note that we also specify a metal shading of 2% and a cell area of $1cm^{2}$. # SolCore calculates the EQE for all three junctions and light-IV showing the relative -# contribution of each sub-cell. We set "kind = 'DA'" to use the depletion -# approximation. +# contribution of each sub-cell. We set "kind = 'PDD'" to use the drift-diffusion +# solver. Note that this will invoke the Fortran PDD solver. To use the Sesame PDD +# solver, we can change the kind to 'sesame_PDD'. We can also mix and match different +# kinds of junction in a single calculation. +# We can also set the surface recombination velocities, where sn +# refers to the surface recombination velocity at the n-type surface, and sp refers +# to the SRV on the p-type side (input units m/s). solar_cell = SolarCell( ARC_layers + @@ -58,7 +62,7 @@ Junction([ Layer(si("25nm"), material=window_material, role='window'), Layer(si("100nm"), material=top_cell_n_material, role='emitter'), - Layer(si("400nm"), material=top_cell_p_material, role='base'), + Layer(si("300nm"), material=top_cell_p_material, role='base'), ], sn=1e4, sp=1e4, kind='PDD'), Junction([Layer(si("200nm"), material=mid_cell_n_material, role='emitter'), Layer(si("3000nm"), material=mid_cell_p_material, role='base'), @@ -70,43 +74,51 @@ shading=0.02, cell_area=1 * 1 / 1e4) # Choose wavelength range (in m): -wl = np.linspace(280, 1850, 200) * 1e-9 +wl = np.linspace(280, 1850, 100) * 1e-9 # Calculate the EQE for the solar cell: -# solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl, -# # 'da_mode': 'bvp', -# 'optics_method': 'TMM' -# }) -# # we pass options to use the TMM optical method to calculate realistic R, A and T -# # values with interference in the ARC (and semiconductor) layers. We can also choose -# # which solver mode to use for the depletion approximation. The default is 'green', -# # which uses the (faster) Green's function method. The other method is 'bvp'. -# -# plt.figure(1) -# plt.plot(wl * 1e9, solar_cell[0 + len(ARC_layers)].eqe(wl) * 100, 'b', label='GaInP QE') -# plt.plot(wl * 1e9, solar_cell[1 + len(ARC_layers)].eqe(wl) * 100, 'g', label='GaAs QE') -# plt.plot(wl * 1e9, solar_cell[2 + len(ARC_layers)].eqe(wl) * 100, 'r', label='Ge QE') -# plt.plot(wl * 1e9, solar_cell[0 + len(ARC_layers)].layer_absorption * 100, '--b', -# label='GaInP Abs.') -# plt.plot(wl * 1e9, solar_cell[1 + len(ARC_layers)].layer_absorption * 100, '--g', -# label='GaAs Abs.') -# plt.plot(wl * 1e9, solar_cell[2 + len(ARC_layers)].layer_absorption * 100, '--r', -# label='Ge Abs.') -# -# plt.plot(wl*1e9, solar_cell.reflected*100, '--k', label="Reflected") -# plt.legend() -# plt.ylim(0, 100) -# plt.ylabel('EQE (%)') -# plt.xlabel('Wavelength (nm)') -# plt.show() +solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl, + 'optics_method': 'TMM' + }) + +# we pass options to use the TMM optical method to calculate realistic R, A and T +# values with interference in the ARC (and semiconductor) layers. +# Plot the EQE and absorption of the individual junctions. Note that we can access +# the properties of the first junction (ignoring any other layers, such as the ARC, +# which are not part of a junction) using solar_cell(0), and the second junction using +# solar_cell(1), etc. + +plt.figure(1) +plt.plot(wl * 1e9, solar_cell(0).eqe(wl) * 100, 'b', label='GaInP QE') +plt.plot(wl * 1e9, solar_cell(1).eqe(wl) * 100, 'g', label='GaAs QE') +plt.plot(wl * 1e9, solar_cell(2).eqe(wl) * 100, 'r', label='Ge QE') +plt.fill_between(wl * 1e9, solar_cell(0).layer_absorption * 100, 0, alpha=0.3, + label='GaInP Abs.', color='b') +plt.fill_between(wl * 1e9, solar_cell(1).layer_absorption * 100, 0, alpha=0.3, + label='GaAs Abs.', color='g') +plt.fill_between(wl * 1e9, solar_cell(2).layer_absorption * 100, 0, alpha=0.3, + label='Ge Abs.', color='r') + +plt.plot(wl*1e9, 100*(1-solar_cell.reflected), '--k', label="100 - Reflectectivity") +plt.legend() +plt.ylim(0, 100) +plt.ylabel('EQE (%)') +plt.xlabel('Wavelength (nm)') +plt.show() # Set up the AM0 (space) solar spectrum am0 = LightSource(source_type='standard',version='AM0',x=wl, output_units='photon_flux_per_m') +# Set up the voltage range for the overall cell (at which the total I-V will be +# calculated) as well as the internal voltages which are used to calculate the results +# for the individual junctions. The range of the internal_voltages should generally +# be wider than that for the voltages. + # this is an n-p cell, so we need to scan negative voltages -V = np.linspace(-3, 0, 50) -internal_voltages = np.linspace(-4, 1.7, 100) + +V = np.linspace(-2.8, 0, 50) +internal_voltages = np.linspace(-3, 1.7, 120) # Calculate the current-voltage relationship under illumination: @@ -118,22 +130,23 @@ 'optics_method': 'TMM', 'mpp': True, }) + # We pass the same options as for solving the EQE, but also set 'light_iv' and 'mpp' to # True to indicate we want the IV curve under illumination and to find the maximum # power point (MPP). We also pass the AM0 light source and voltages created above. plt.figure(2) -plt.plot(abs(V), -solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell') -plt.plot(abs(V), solar_cell[0 + len(ARC_layers)].iv(V)/10, 'b', label='InGaP sub-cell') -plt.plot(abs(V), solar_cell[1 + len(ARC_layers)].iv(V)/10, 'g', label='GaAs sub-cell') -plt.plot(abs(V), solar_cell[2 + len(ARC_layers)].iv(V)/10, 'r', label='Ge sub-cell') +plt.plot(-solar_cell.iv['IV'][0], -solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell') +plt.plot(-internal_voltages, solar_cell(0).iv(internal_voltages)/10, 'b', label='InGaP sub-cell') +plt.plot(-internal_voltages, solar_cell(1).iv(internal_voltages)/10, 'g', label='GaAs sub-cell') +plt.plot(-internal_voltages, solar_cell(2).iv(internal_voltages)/10, 'r', label='Ge sub-cell') plt.text(0.5,30,f'Jsc= {abs(solar_cell.iv.Isc/10):.2f} mA.cm' + r'$^{-2}$') plt.text(0.5,28,f'Voc= {abs(solar_cell.iv.Voc):.2f} V') plt.text(0.5,26,f'FF= {solar_cell.iv.FF*100:.2f} %') plt.text(0.5,24,f'Eta= {solar_cell.iv.Eta*100:.2f} %') plt.legend() -plt.ylim(0, 35) +plt.ylim(-10, 35) plt.xlim(0, 3) plt.ylabel('Current (mA/cm$^2$)') plt.xlabel('Voltage (V)') diff --git a/examples/MJ_solar_cell_using_DA.py b/examples/MJ_solar_cell_using_DA.py index be983fbd..85c0d060 100755 --- a/examples/MJ_solar_cell_using_DA.py +++ b/examples/MJ_solar_cell_using_DA.py @@ -27,16 +27,6 @@ hole_diffusion_length=si("200nm")) top_cell_p_material = InGaP(In=0.49, Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("2um")) - -# For convenience we will set the carrier transport properties to default values for -# all the materials in one place. To do that we add the materials to a list called -# all_materials: - -all_materials = [] -all_materials.append(window_material) -all_materials.append(top_cell_n_material) -all_materials.append(top_cell_p_material) - # MID CELL - GaAs GaAs = material("GaAs") @@ -45,9 +35,6 @@ mid_cell_p_material = GaAs(In=0.01, Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("5um")) -all_materials.append(mid_cell_n_material) -all_materials.append(mid_cell_p_material) - # BOTTOM CELL - Ge Ge = material("Ge") @@ -56,24 +43,13 @@ bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("50um")) -all_materials.append(bot_cell_n_material) -all_materials.append(bot_cell_p_material) - -# We now set the electronic properties of these materials to some default values. -# Although in practice there will be differences, for this example we set them to -# reasonable generic values for simplicity. - -for mat in all_materials: - mat.hole_mobility = 5e-2 - mat.electron_mobility = 3.4e-3 - mat.hole_mobility = 3.4e-3 - mat.electron_mobility = 5e-2 - # Now that the layers are configured, we can now assemble the triple junction solar # cell. Note that we also specify a metal shading of 2% and a cell area of $1cm^{2}$. # SolCore calculates the EQE for all three junctions and light-IV showing the relative # contribution of each sub-cell. We set "kind = 'DA'" to use the depletion -# approximation. +# approximation. We can also set the surface recombination velocities, where sn +# refers to the surface recombination velocity at the n-type surface, and sp refers +# to the SRV on the p-type side. solar_cell = SolarCell( ARC_layers + @@ -81,13 +57,13 @@ Junction([Layer(si("25nm"), material=window_material, role='window'), Layer(si("100nm"), material=top_cell_n_material, role='emitter'), Layer(si("400nm"), material=top_cell_p_material, role='base'), - ], sn=1, sp=1, kind='DA'), + ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA'), Junction([Layer(si("200nm"), material=mid_cell_n_material, role='emitter'), Layer(si("3000nm"), material=mid_cell_p_material, role='base'), - ], sn=1, sp=1, kind='DA'), + ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA'), Junction([Layer(si("400nm"), material=bot_cell_n_material, role='emitter'), Layer(si("100um"), material=bot_cell_p_material, role='base'), - ], sn=1, sp=1, kind='DA') + ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA') ], shading=0.02, cell_area=1 * 1 / 1e4) @@ -96,7 +72,7 @@ # Calculate the EQE for the solar cell: solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl, - # 'da_mode': 'bvp', + 'da_mode': 'green', 'optics_method': 'TMM' }) # we pass options to use the TMM optical method to calculate realistic R, A and T @@ -104,28 +80,41 @@ # which solver mode to use for the depletion approximation. The default is 'green', # which uses the (faster) Green's function method. The other method is 'bvp'. +# Plot the EQE and absorption of the individual junctions. Note that we can access +# the properties of the first junction (ignoring any other layers, such as the ARC, +# which are not part of a junction) using solar_cell(0), and the second junction using +# solar_cell(1), etc. + plt.figure(1) -plt.plot(wl * 1e9, solar_cell[0 + len(ARC_layers)].eqe(wl) * 100, 'b', label='GaInP QE') -plt.plot(wl * 1e9, solar_cell[1 + len(ARC_layers)].eqe(wl) * 100, 'g', label='GaAs QE') -plt.plot(wl * 1e9, solar_cell[2 + len(ARC_layers)].eqe(wl) * 100, 'r', label='Ge QE') -plt.plot(wl * 1e9, solar_cell[0 + len(ARC_layers)].layer_absorption * 100, '--b', - label='GaInP Abs.') -plt.plot(wl * 1e9, solar_cell[1 + len(ARC_layers)].layer_absorption * 100, '--g', - label='GaAs Abs.') -plt.plot(wl * 1e9, solar_cell[2 + len(ARC_layers)].layer_absorption * 100, '--r', - label='Ge Abs.') - -plt.plot(wl*1e9, solar_cell.reflected*100, '--k', label="Reflected") +plt.plot(wl * 1e9, solar_cell(0).eqe(wl) * 100, 'b', label='GaInP QE') +plt.plot(wl * 1e9, solar_cell(1).eqe(wl) * 100, 'g', label='GaAs QE') +plt.plot(wl * 1e9, solar_cell(2).eqe(wl) * 100, 'r', label='Ge QE') +plt.fill_between(wl * 1e9, solar_cell(0).layer_absorption * 100, 0, alpha=0.3, + label='GaInP Abs.', color='b') +plt.fill_between(wl * 1e9, solar_cell(1).layer_absorption * 100, 0, alpha=0.3, + label='GaAs Abs.', color='g') +plt.fill_between(wl * 1e9, solar_cell(2).layer_absorption * 100, 0, alpha=0.3, + label='Ge Abs.', color='r') + +plt.plot(wl*1e9, 100*(1-solar_cell.reflected), '--k', label="100 - Reflectivity") plt.legend() plt.ylim(0, 100) plt.ylabel('EQE (%)') plt.xlabel('Wavelength (nm)') plt.show() -# Set up the AM0 (space) solar spectrum +# Set up the AM0 (space) solar spectrum for the light I-V calculation: am0 = LightSource(source_type='standard',version='AM0',x=wl, output_units='photon_flux_per_m') + +# Set up the voltage range for the overall cell (at which the total I-V will be +# calculated) as well as the internal voltages which are used to calculate the results +# for the individual junctions. The range of the internal_voltages should generally +# be wider than that for the voltages. + +# this is an n-p cell, so we need to scan negative voltages + V = np.linspace(-3, 0, 300) internal_voltages = np.linspace(-4, 2, 400) @@ -138,19 +127,19 @@ 'wavelength': wl, 'optics_method': 'TMM', 'mpp': True, - # 'da_mode': 'bvp' }) + # We pass the same options as for solving the EQE, but also set 'light_iv' and 'mpp' to # True to indicate we want the IV curve under illumination and to find the maximum # power point (MPP). We also pass the AM0 light source and voltages created above. plt.figure(2) plt.plot(abs(V), -solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell') -plt.plot(abs(V), solar_cell[0 + len(ARC_layers)].iv(V)/10, 'b', label='InGaP sub-cell') -plt.plot(abs(V), solar_cell[1 + len(ARC_layers)].iv(V)/10, 'g', label='GaAs sub-cell') -plt.plot(abs(V), solar_cell[2 + len(ARC_layers)].iv(V)/10, 'r', label='Ge sub-cell') -plt.text(0.5,30,f'Jsc= {solar_cell.iv.Isc/10:.2f} mA.cm' + r'$^{-2}$') -plt.text(0.5,28,f'Voc= {solar_cell.iv.Voc:.2f} V') +plt.plot(abs(V), solar_cell(0).iv(V)/10, 'b', label='InGaP sub-cell') +plt.plot(abs(V), solar_cell(1).iv(V)/10, 'g', label='GaAs sub-cell') +plt.plot(abs(V), solar_cell(2).iv(V)/10, 'r', label='Ge sub-cell') +plt.text(0.5,30,f'Jsc= {abs(solar_cell.iv.Isc/10):.2f} mA.cm' + r'$^{-2}$') +plt.text(0.5,28,f'Voc= {abs(solar_cell.iv.Voc):.2f} V') plt.text(0.5,26,f'FF= {solar_cell.iv.FF*100:.2f} %') plt.text(0.5,24,f'Eta= {solar_cell.iv.Eta*100:.2f} %') diff --git a/examples/Si_cell.py b/examples/Si_cell.py index c4fe2d28..b43663f6 100644 --- a/examples/Si_cell.py +++ b/examples/Si_cell.py @@ -2,23 +2,21 @@ from solcore.state import State from solcore.solar_cell_solver import solar_cell_solver import numpy as np -import os from solcore.light_source import LightSource -from solcore.absorption_calculator import search_db from scipy.special import erfc + import matplotlib.pyplot as plt from solcore import material, si -d_bulk = 130e-6 +d_bulk = 100e-6 Air = material("Air")() -MgF2_pageid = str(search_db(os.path.join("MgF2", "Rodriguez-de Marcos"))[0][0]) -MgF2 = material(MgF2_pageid, nk_db=True)() -Ag = material("Ag_Jiang")() +MgF2 = material("MgF2")() +Ag = material("Ag")() Si_pn = material("Si")(electron_mobility=si("1e4cm2"), hole_mobility=si("1e3cm2"), - electron_minority_lifetime=0.01, hole_minority_lifetime=0.01) + electron_minority_lifetime=0.001, hole_minority_lifetime=0.001) wavelengths = np.linspace(280, 1200, 100)*1e-9 @@ -36,13 +34,13 @@ options.T = 298 options.light_source = LightSource(source_type="standard", - version="AM1.5g", x=options.wavelength, output_units="photon_flux_per_m") + version="AM1.5g", x=options.wavelength) -options.voltages = np.linspace(0, 0.8, 30) +options.voltages = np.linspace(0, 0.8, 40) options.internal_voltages = options.voltages options.mpp = True -options.recalculate_absorption = True options.no_back_reflection = False +options.position = 10e-9 nD = si("1e20cm-3") nA = si("1e19cm-3") @@ -59,12 +57,10 @@ def doping_profile_func(x): return doping_profile + doping_profile_rear + bulk_doping -L_cm = d_bulk*100 # length of the system in the x-direction [cm] - -surface_recomb = dict(sn_front=si('5e3 cm s-1'), # important - sp_front=si('1e4 cm s-1'), - sn_rear=si('1e4 cm s-1'), - sp_rear=si('1e4 cm s-1')) # important +surface_recomb = dict(sp_e=si('100 cm s-1'), # important (minority carrier at front) + sp_h=si('100 cm s-1'), + sn_e=si('100 cm s-1'), + sn_h=si('100 cm s-1')) # important (minority carrier at rear) Si_junction = [Junction([Layer(d_bulk, Si_pn)], doping_profile=doping_profile_func, kind='sesame_PDD', @@ -72,7 +68,6 @@ def doping_profile_func(x): ) ] - Si_cell = SolarCell(front_materials + Si_junction + back_materials, @@ -90,7 +85,7 @@ def doping_profile_func(x): ax.stackplot(wavelengths * 1e9, 100 * result_stack[::-1], linewidth=0.5, alpha=0.5, labels=['MgF2 (rear)', 'TCO (rear)', 'Si bulk', 'TCO (front)', 'MgF2 (front)', 'Reflection']) -ax.plot(wavelengths * 1e9, 100 * Si_cell(0).eqe(wavelengths), '-r', linewidth=2, +ax.plot(wavelengths * 1e9, 100 * Si_cell(0).eqe(wavelengths), '-k', linewidth=2, label='EQE') ax.set_xlim(280, 1200) @@ -135,7 +130,7 @@ def doping_profile_func(x): ax2.text(0.1, 0.9 * jsc, r'= {:.2f} mA/cm$^2$'.format(jsc)) ax2.text(0.1, 0.8 * jsc, r'= {:.3f} V'.format(Si_cell.iv.Voc)) ax2.text(0.1, 0.7 * jsc, '= {:.2f} %'.format(Si_cell.iv.FF * 100)) -ax2.text(0.1, 0.6 * jsc, r'= {:.2f} %'.format(Si_cell.iv.Pmpp / 10)) +ax2.text(0.1, 0.6 * jsc, r'= {:.2f} %'.format(Si_cell.iv.Eta * 100)) ax2.text(0.1, 0.5 * jsc, r'= {:.2f} mA/cm$^2$'.format( Si_cell.iv.Impp / 10)) ax2.text(0.1, 0.4 * jsc, r'= {:.3f} V'.format(Si_cell.iv.Vmpp)) ax2.grid(which='major', alpha=0.35) @@ -144,3 +139,59 @@ def doping_profile_func(x): plt.tight_layout() plt.show() + +# Scan through lifetimes + +lifetime_exp = np.arange(-4, 0, 1) + +lifetimes = 10.0**lifetime_exp + +cell_results = np.zeros(([len(lifetimes), 4])) # save efficiency, FF, Voc, Jsc + +for i1, lt in enumerate(lifetimes): + + options.recalculate_absorption = True + + Si_pn = material("Si")(electron_mobility=si("1e4cm2"), hole_mobility=si("1e3cm2"), + electron_minority_lifetime=lt, hole_minority_lifetime=lt) + + Si_junction = [Junction([Layer(d_bulk, Si_pn)], + doping_profile=doping_profile_func, kind='sesame_PDD', + **surface_recomb)] + + Si_cell = SolarCell(front_materials + + Si_junction + + back_materials, + shading=0.02, + substrate=Ag, + ) + + solar_cell_solver(Si_cell, 'iv', options) + + cell_results[i1] = np.array([100*Si_cell.iv.Eta, 100*Si_cell.iv.FF, Si_cell.iv.Voc, Si_cell.iv.Isc/10]) + + print(lt, Si_cell.iv.Eta, Si_cell.iv.FF) + + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.5)) + +ax1.plot(cell_results[:, 2], cell_results[:, 1], 'ko') + +for i, lt in enumerate(lifetimes): + ax1.annotate(str(lt), (cell_results[i, 2] - 0.001, cell_results[i, 1]), ha='right') + +ax1.set_xlabel(r'V$_{oc}$ (V)') +ax1.set_ylabel('FF (%)') +ax1.set_xlim(0.98*np.min(cell_results[:, 2]), 1.01*np.max(cell_results[:, 2])) + +ax2.semilogx(lifetimes, cell_results[:, 0], 'o', color='k') +ax2.set_ylabel('Efficiency (%)') +ax3 = ax2.twinx() +ax3.plot(lifetimes, cell_results[:, 3], 'o', color='r', markerfacecolor='none') +ax3.set_ylim(32, 35) +ax3.set_ylabel(r'$J_{sc}$ (mA/cm$^2$)', color='r') + +ax2.set_xlabel(r'$\tau$ (s)') + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/examples/sesame_PDD_example.py b/examples/sesame_PDD_example.py deleted file mode 100644 index d80b3dce..00000000 --- a/examples/sesame_PDD_example.py +++ /dev/null @@ -1,305 +0,0 @@ -from solcore import material, si -from solcore.solar_cell import SolarCell, Junction, Layer -from solcore.state import State -from solcore.solar_cell_solver import solar_cell_solver -from solcore.light_source import LightSource -import numpy as np -import matplotlib.pyplot as plt - -options = State() -# options.wavelength = np.linspace(280, 700, 20)*1e-9 -# options.optics_method = 'TMM' -# options.light_iv = True -# options.light_source = LightSource(source_type="standard", -# version="AM1.5g", x=options.wavelength, output_units="photon_flux_per_m") -# options.voltages = np.linspace(0, 2, 40) -# # options.minimum_spacing = 1e-09 -# # options.maximum_spacing = 1e-8 -# -# T = 293 -# -# add_args = {'relative_permittivity': 10, 'electron_minority_lifetime': 5e-6, -# 'hole_minority_lifetime': 5e-6, -# 'electron_auger_recombination': 1e-45, -# 'hole_auger_recombination': 1e-45} -# -# ARC = material('Si3N4')() -# window = material('AlGaAs')(T=T, Na=5e24, Al=0.8, **add_args) -# p_AlGaAs = material('AlGaAs')(T=T, Na=1e24, Al=0.4, **add_args) -# n_AlGaAs = material('AlGaAs')(T=T, Nd=8e22, Al=0.4, **add_args) -# bsf = material('AlGaAs')(T=T, Nd=2e24, Al=0.6, **add_args) -# -# junction = Junction([Layer(width=si('30nm'), material=window, role="Window"), -# Layer(width=si('1000nm'), material=p_AlGaAs, role="Emitter"), -# Layer(width=si('150nm'), material=n_AlGaAs, role="Base"), -# Layer(width=si('200nm'), material=bsf, role="BSF")], sn=1e6, sp=1e6, T=T, kind='sesame_PDD') -# -# widths = [layer.width for layer in junction] -# -# # x = np.concatenate((np.linspace(0,3e-7, 400, endpoint=False), -# # np.linspace(3e-7, np.sum(widths) - 100e-9, 100, endpoint=False), -# # np.linspace(np.sum(widths) - 100e-9, np.sum(widths), 200))) -# # junction.mesh = x -# -# solar_cell = SolarCell( -# [Layer(60e-0, ARC), junction] -# ) -# -# solar_cell_solver(solar_cell, 'iv', options) -# solar_cell_solver(solar_cell, 'qe', options) -# -# # # qe_sesame(solar_cell[1], options) -# # -# # iv_sesame(solar_cell[1], options) -# -# junction_old = Junction([Layer(width=si('30nm'), material=window, role="Window"), -# Layer(width=si('150nm'), material=p_AlGaAs, role="Emitter"), -# Layer(width=si('1000nm'), material=n_AlGaAs, role="Base"), -# Layer(width=si('200nm'), material=bsf, role="BSF")], sn=1e6, sp=1e6, T=T, -# kind='PDD') -# -# -# solar_cell_old = SolarCell( -# [Layer(60e-0, ARC), junction_old] -# ) -# -# # solar_cell_solver(solar_cell_old, 'iv', options) -# # -# plt.figure() -# plt.plot(options.voltages, solar_cell[1].iv(options.voltages), 'o-', label='IV') -# # plt.plot(options.voltages, solar_cell_old[1].iv(options.voltages), 'o-', label='IV') -# # plt.xlim(-2, 1.8) -# plt.ylim(-200, 200) -# plt.show() -# -# plt.figure() -# plt.plot(options.wavelength*1e9, solar_cell.absorbed) -# # plt.plot(options.wavelength*1e9, solar_cell_old.absorbed) -# plt.show() -# -# plt.figure() -# plt.plot(options.wavelength*1e9, solar_cell[1].eqe(options.wavelength)) -# plt.plot(options.wavelength*1e9, solar_cell[1].iqe(options.wavelength)) -# plt.show() -# -# solar_cell[1].absorbed(solar_cell[1].mesh) -# -# plt.figure() -# plt.plot(junction.mesh, 'x') -# plt.show() -# -# - -##### - -import os -from solcore.constants import q -from solcore.light_source import LightSource -from solcore.interpolate import interp1d -from solcore.absorption_calculator import search_db -from scipy.special import erfc -import matplotlib.pyplot as plt -from rayflare.textures import regular_pyramids - -from rayflare.transfer_matrix_method import tmm_structure - -from rayflare.textures import regular_pyramids, rough_pyramids -from solcore import material, si -from rayflare.options import default_options -from solcore.structure import Layer -from rayflare.ray_tracing import rt_structure -from rayflare.utilities import make_absorption_function -from os.path import exists - -# Note Sesame assumes all quantities of length are input in units of cm. Other assumed -# input units are: time in s, energy in eV. - -# LONGi paper: -# To further increase JSC, a 150-nm-thick MgF2 film was evaporated on the front TCO layer -# as a second antireflective coating. For 26.81% cell, an additional 120-nm-thick MgF2/150-nm-thick -# Ag stack was evaporated on the rear TCO layer, which means this cell is a monofacial solar cell. - -options = State() -options.wavelength = np.linspace(280, 1200, 20)*1e-9 -options.optics_method = 'TMM' -options.light_iv = True -options.light_source = LightSource(source_type="standard", - version="AM1.5g", x=options.wavelength, output_units="photon_flux_per_m") -options.voltages = np.linspace(0, 0.8, 20) -options.mpp = True - -# ray-tracing for optics first -d_bulk = 130e-6 - -# ITO_front = material("ITO_front_transparent")() -ITO_front = material("ITO_front_transparent")() -ITO_back = material("ITO_back")() -Ag = material("Ag_Jiang")() -aSi_i = material("aSi_i")() -aSi_p = material("aSi_p")() -aSi_n = material("aSi_n")() -Air = material("Air")() -MgF2_pageid = str(search_db(os.path.join("MgF2", "Rodriguez-de Marcos"))[0][0]) -MgF2 = material(MgF2_pageid, nk_db=True)() -Ag = material("Ag_Jiang")() - -Si_pn = material("Si")(Nc=3e25, Nv=1e25, electron_mobility=si("1e4cm2"), hole_mobility=si("1e3cm2"), - electron_minority_lifetime=0.05, hole_minority_lifetime=0.05, - radiative_recombination=si("1.89e-15cm3"), - electron_auger_recombination=si("3e-31cm6"), hole_auger_recombination=si("1e-31cm6")) - -front_materials = [Layer(110e-9, MgF2), Layer(40 * 1e-9, ITO_front), - Layer(1e-9, aSi_i), Layer(1e-9, aSi_n)] -back_materials = [Layer(2e-9, aSi_i), Layer(2e-9, aSi_p), Layer(100e-9, ITO_back), - Layer(120e-9, MgF2), Layer(150e-9, Ag)] - -front_width = np.sum([layer.width for layer in front_materials]) -back_width = np.sum([layer.width for layer in back_materials]) - -L = 130e-6 - -# rear junction - -nD = si("1e20cm-3") -nA = si("1e19cm-3") -bulk_doping = si("1e15cm-3") # n type bulk -def doping_profile_func(x): - - L = d_bulk - - doping_profile = - nA * erfc(x/150e-9) # characteristic depth of 150 nm - - doping_profile_rear = nD * erfc((L - x)/200e-9) # characteristic depth of 200 nm - - return doping_profile + doping_profile_rear + bulk_doping - -L_cm = d_bulk*100 # length of the system in the x-direction [cm] - -# Mesh -x = np.concatenate((np.linspace(0,2e-4, 2000, endpoint=False), - np.linspace(2e-4, L_cm - 2e-4, 2000, endpoint=False), - np.linspace(L_cm - 2e-4, L_cm, 2000))) - -options.position = np.concatenate((np.linspace(0, front_width, 100), front_width + x/100, - np.linspace(front_width + d_bulk, - front_width + d_bulk + back_width, 100))) - -options.minimum_spacing = 5e-9 - -Si_junction = [Junction([Layer(d_bulk, Si_pn)], - doping_profile=doping_profile_func, kind='sesame_PDD', - # mesh=x/100 - )] - -Si_cell = SolarCell(front_materials + - Si_junction + - back_materials) - -Si_ind = len(front_materials) - -solar_cell_solver(Si_cell, 'qe', options) - -plt.figure() -plt.plot(options.wavelength*1e9, Si_cell[Si_ind].layer_absorption, '--', label='Si absorption') -plt.plot(options.wavelength*1e9, Si_cell[Si_ind].eqe(options.wavelength), '-k', label='Si EQE') -plt.plot(options.wavelength*1e9, Si_cell[Si_ind].iqe(options.wavelength), '-r', label='Si IQE') -plt.legend() -plt.show() - -solar_cell_solver(Si_cell, 'iv', options) - -jsc = Si_cell.iv.Isc/10 - -plt.figure() -plt.plot(options.voltages, -Si_cell[Si_ind].iv(options.voltages)/10, 'o-', label='IV') -# plt.plot(options.voltages, solar_cell_old[1].iv(options.voltages), 'o-', label='IV') -# plt.xlim(-2, 1.8) -plt.text(0.02, 0.9*jsc, r'Jsc = {:.2f} mA/cm$^2$'.format(jsc)) -plt.text(0.02, 0.8*jsc, 'Voc = {:.3f} V'.format(Si_cell.iv.Voc)) -plt.text(0.02, 0.7*jsc, 'FF = {:.3f} %'.format(Si_cell.iv.FF*100)) -plt.text(0.02, 0.6*jsc, 'Eff = {:.3f} %'.format(Si_cell.iv.Eta*100)) -plt.text(0.02, 0.5*jsc, r'Jmpp = {:.3f} mA/cm$^2$'.format(Si_cell.iv.Impp/10)) -plt.text(0.02, 0.4*jsc, 'Vmpp = {:.3f} V'.format(Si_cell.iv.Vmpp)) -plt.ylim(0, 1.1*jsc) -plt.xlim(np.min(options.voltages), np.max(options.voltages)) -plt.xlabel('Voltage (V)') -plt.ylabel('Current density (mA/cm$^2$)') -plt.show() - -# plt.figure() -# plt.plot(x*1e4, Si_cell[Si_ind].sesame_sys.g, 'o') -# plt.xlim(0, 0.1) -# plt.show() -# -# plt.figure() -# plt.semilogy(Si_junction[0].mesh*1e6, Si_cell[Si_ind].absorbed(Si_junction[0].mesh), 'o', markerfacecolor='none') -# plt.xlim(0, 1) -# plt.ylim(1e-3,1e8) -# plt.show() - -# -# x = np.linspace(0, d_bulk, int(1e7)) -# int_A = np.trapz(Si_cell[Si_ind].absorbed(x), x, axis=0) -# def doping_profile_test(x): -# -# L = 8000e-9 -# -# doping_profile = - 1e20 * erfc(x / 150e-9) # characteristic depth of 150 nm -# -# doping_profile_2 = - 5e19 * erfc(x/400e-9) # characteristic depth of 150 nm -# # doping_profile = 0 -# # doping_profile_2 = 0 -# doping_profile_rear = 5e19 * erfc(((L/2)- x)/200e-9) # characteristic depth of 200 nm -# -# doping_profile_3 = 1e20 * erfc((L - x) / 150e-9) # characteristic depth of 150 nm -# -# return (doping_profile + doping_profile_2 + doping_profile_3 + -# doping_profile_rear + 1e16) -# -# x_test = np.linspace(0, 8000, 1000)*1e-9 -# -# doping_test = doping_profile_test(x_test) -# -# diff_rho = np.abs(np.gradient(doping_test, x_test)) -# -# -# -# -# doping_change = np.where(diff_rho/np.max(diff_rho) > 1e-3)[0] -# -# # identify continuous stretches of high doping change, and where there are gaps: -# -# change_point = np.where(np.diff(doping_change) != 1)[0] -# -# last_points = doping_change[change_point] -# first_points = doping_change[change_point + 1] -# -# first_points = np.insert(first_points, 0, doping_change[0]) -# last_points = np.append(last_points, doping_change[-1]) -# -# first_x = x_test[first_points] -# last_x = x_test[last_points] -# -# plt.figure() -# plt.plot(x_test*1e9, doping_test/np.max(doping_test), 'o-') -# plt.plot(x_test*1e9, diff_rho/np.max(diff_rho), 'o--') -# plt.plot(x_test[first_points]*1e9, [0.5]*len(first_points), 'xr') -# plt.plot(x_test[last_points]*1e9, [0.5]*len(last_points), 'xk') -# plt.show() -# -# current_mesh = x_test -# -# minimum_spacing = 1e-9 -# -# for i1 in range(len(first_points)): -# # want to make mesh denser where there are large changes in doping -# -# mesh_above = current_mesh[current_mesh < first_x[i1]] -# mesh_below = current_mesh[current_mesh >= last_x[i1]] -# mesh_dense = np.arange(first_x[i1], last_x[i1] - 1e-12, minimum_spacing) # last point not included -# -# current_mesh = np.unique(np.concatenate((mesh_above, mesh_dense, mesh_below))) -# -# plt.figure() -# plt.plot(current_mesh*1e9) -# plt.show() \ No newline at end of file diff --git a/examples/sesame_testing.py b/examples/sesame_testing.py deleted file mode 100644 index 037bfd81..00000000 --- a/examples/sesame_testing.py +++ /dev/null @@ -1,478 +0,0 @@ -import sesame -import numpy as np -import os -from solcore.constants import q -from solcore.light_source import LightSource -from solcore.interpolate import interp1d -from solcore.absorption_calculator import search_db -from scipy.special import erfc -import matplotlib.pyplot as plt -from rayflare.textures import regular_pyramids - -from rayflare.transfer_matrix_method import tmm_structure - -from rayflare.textures import regular_pyramids, rough_pyramids -from solcore import material, si -from rayflare.options import default_options -from solcore.structure import Layer -from rayflare.ray_tracing import rt_structure -from rayflare.utilities import make_absorption_function -from os.path import exists - -# Note Sesame assumes all quantities of length are input in units of cm. Other assumed -# input units are: time in s, energy in eV. - -# LONGi paper: -# To further increase JSC, a 150-nm-thick MgF2 film was evaporated on the front TCO layer -# as a second antireflective coating. For 26.81% cell, an additional 120-nm-thick MgF2/150-nm-thick -# Ag stack was evaporated on the rear TCO layer, which means this cell is a monofacial solar cell. - -force_recalculate = True - -step_size = 20 - -shading = 0.02 - -wls = np.arange(280, 1200, step_size) * 1e-9 - -wls_ls = np.linspace(250, 4500, 4000) -base_light = LightSource(source_type="standard", - version="AM1.5g", x=wls_ls, output_units="photon_flux_per_nm") - -base_flux = base_light.spectrum()[1] - -# ray-tracing for optics first -d_bulk = 130e-6 - -# photons at wavelength shorter than the singlet energy S1 are assumed to be absorbed, -# and turned into two photons through singlet fission. Assuming thermally neutral SF, these -# two photons will have half the energy/two times the wavelength of S1. - -absorption_data = f'rt_data_{step_size}3.npy' - -# ITO_front = material("ITO_front_transparent")() -ITO_front = material("ITO_front_transparent")() -ITO_back = material("ITO_back")() -Ag = material("Ag_Jiang")() -aSi_i = material("aSi_i")() -aSi_p = material("aSi_p")() -aSi_n = material("aSi_n")() -Si = material("Si")() -Air = material("Air")() -SiNx = material("SiNx_Ox")() -MgF2_pageid = str(search_db(os.path.join("MgF2", "Rodriguez-de Marcos"))[0][0]) -MgF2 = material(MgF2_pageid, nk_db=True)() -Ag = material("Ag_Jiang")() - - -plt.figure() -plt.plot(wls*1e9, SiNx.n(wls)) -plt.plot(wls*1e9, SiNx.k(wls)) -plt.plot(wls*1e9, ITO_front.n(wls), '--') -plt.plot(wls*1e9, ITO_front.k(wls), '--') -plt.plot(wls*1e9, ITO_back.n(wls), '-.') -plt.plot(wls*1e9, ITO_back.k(wls), '-.') -plt.show() -# build structure - -options = default_options() -options.wavelengths = wls -options.pol = "u" -options.I_thresh = 1e-3 -options.randomize_surface = True -options.n_jobs = -3 -options.periodic = True -options.project_name = "basic_cell_model" -options.nx = 20 -options.ny = options.nx -options.n_rays = 4*options.nx ** 2 -options.lookuptable_angles = 500 -options.depth_spacing = 1e-9 -options.depth_spacing_bulk = 1e-9 -options.theta_in = 0*np.pi/180 -options.phi_in = 0*np.pi/180 -options.parallel = True - -options.theta_in = 0*np.pi/180 - - -front_materials = [Layer(110e-9, MgF2), Layer(40 * 1e-9, ITO_front), - Layer(1e-9, aSi_i), Layer(1e-9, aSi_n)] -back_materials = [Layer(2e-9, aSi_i), Layer(2e-9, aSi_p), Layer(100e-9, ITO_back), - Layer(120e-9, MgF2), Layer(150e-9, Ag)] - -# thickness of front: -d_front = np.sum([x.width*1e9 for x in front_materials]) - -tmm_str = tmm_structure(front_materials + [Layer(d_bulk, Si)] + back_materials, incidence=Air, transmission=Air) - -options.coherent = False -options.coherency_list = ['c']*len(front_materials) + ['i']*(len(back_materials) + 1) - -tmm_result = tmm_str.calculate_profile(options) - -bulk_profile = tmm_result['profile'] - -Si_profile = bulk_profile[:, int(d_front):(int(1e9*d_bulk)+int(d_front)+1)] - - -# plot integrated absorption - this is what we actually need to compare the EQE to see where losses are -# coming from. - -bulk_positions_cm = np.arange(0, d_bulk * 1e2, options.depth_spacing_bulk * 1e2) - -# convert bulk_profile from nm-1 to cm-1 - -profile_func = interp1d(bulk_positions_cm, 1e7 * Si_profile, kind='linear', bounds_error=False, fill_value=0) - -L = d_bulk*100 # length of the system in the x-direction [cm] - -# Mesh -x = np.concatenate((np.linspace(0,2e-4, 2000, endpoint=False), - np.linspace(2e-4, L - 2e-4, 2000, endpoint=False), - np.linspace(L - 2e-4, L, 2000))) - -A_intgr = np.trapz(profile_func(x), x, axis=1) - -Si = material("Si")() - -nD = 1e19 # [cm^-3] -nA = 1e20 # [cm^-3] -bulk_doping = 1e15 # p-type - -doping_profile = nD * erfc(x/150e-7) # characteristic depth of 150 nm - -doping_profile_rear = - nA * erfc((L - x)/200e-7) # characteristic depth of 200 nm - -overall_profile = doping_profile + doping_profile_rear - bulk_doping -# Create a system -sys = sesame.Builder(x) - -# https://www.pveducation.org/es/fotovoltaica/materials/general-properties-of-silicon - -# Dictionary with the material parameters -material_Si = {'Nc': 3e19, 'Nv': 1e19, 'Eg': 1.12, 'affinity': 4.05, 'epsilon': 11.7, - 'mu_e': 1000, 'mu_h': 10000, 'tau_e': 0.05, 'tau_h': 0.05, 'Et': 0, - 'B': 1.89e-15, - 'Cn': 3e-31, 'Cp': 1e-31 - } - -# Auger: https://ui.adsabs.harvard.edu/abs/1977ApPhL..31..346D/abstract#:~:text=The%20Auger%20coefficients%20at%20300,31%20cm6%20s%2D1. - - -# Et: energy levels of the bulk recombination centers [eV] -# B: radiation recombination constant [cm\ :sup:`3`/s], Cn (Cp): Auger - # recombination constant for electrons (holes) [cm\ :sup:`6`/s], - # mass_e (mass_h): effective mass of electrons (holes). All parameters - # can be scalars or callable functions similar to the location - # argument. -# can all be location dependent - - - -# Add the material to the system -sys.add_material(material_Si) - -# fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True) -# ax1.plot(x, overall_profile, label='erfc') -# # ax1.plot(x, sys.rho*sys.scaling.density, label='flat') -# ax1.set_xlim(0, 400e-7) -# ax1.legend() -# -# ax2.plot(x, overall_profile) -# # ax2.plot(x, sys.rho*sys.scaling.density) -# ax2.set_xlim(L - 600e-7, L) -# -# plt.show() - -sys.rho = overall_profile/sys.scaling.density - -# Define Ohmic contacts -sys.contact_type('Ohmic', 'Ohmic') - -# Define the surface recombination velocities for electrons and holes [cm/s] -Sn_left, Sp_left, Sn_right, Sp_right = 1, 1e3, 1e3, 1 -# lower than 1e-6: diverges -sys.contact_S(Sn_left, Sp_left, Sn_right, Sp_right) - -def gfcn_fun_rt(wl_index, flux): - - # profile_func is in m-1, but as a function of cm - - def gcfn(x, y): - return flux*profile_func(x)[wl_index] - - return gcfn - -test_rt = gfcn_fun_rt(8, 1e20) -rt_data = test_rt(x, 0) - -eqe = np.zeros_like(wls) - -for i1, wl in enumerate(wls): - - flux = 1e20 - print(wl) - sys.generation(gfcn_fun_rt(i1, flux)) - print(gfcn_fun_rt(i1, flux)(0, 0)) - print(gfcn_fun_rt(i1, flux)(1e-6, 0)) - voltages = [0] - - if i1 == 0: - guess = None - - else: - guess = {key: result[key][0,:] for key in result.keys()} - - j, result = sesame.IVcurve(sys, voltages, guess=guess) - - # A[i1] = 1 - np.exp(-alpha*L) - # eqe[i1] = j/(q*flux) - - eqe[i1] = j/(q*flux) - - -A = tmm_result['A_per_layer'][:, len(front_materials)] -iqe = eqe/A - -# convert dimensionless current to dimension-ful current -eqe = eqe * sys.scaling.current -iqe = iqe * sys.scaling.current - -plt.figure() -plt.plot(wls*1e9, eqe,'-b', label='EQE') -plt.plot(wls*1e9, iqe,'--r', label='IQE') -plt.plot(wls*1e9, A_intgr, '--k', label='A') -plt.xlabel('Wavelength (nm)') -plt.ylabel('EQE') -plt.legend() -plt.grid() -plt.show() - -# def gen_wl(x): -# return 1e-2*Si.alpha(wls)[None,:]*1e-2 * np.exp(-1e-2*Si.alpha(wls)[None,:]* x[:, -# None]) - -gen_wl = profile_func - -# gg = light_source.spectrum()[1][None,:] * gen_wl(x) -# photon flux m-2 / m, generation cm-1 -gg = base_light.spectrum(wls, output_units="photon_flux_per_m")[1][:, None] * gen_wl(x) -# g_vs_z = np.trapz(gg, wls, axis=1) / 1e4 # m-2 cm-1 -> cm-3 -g_vs_z = (1-shading)*np.trapz(gg, wls, axis=0) / 1e4 -g_vs_z[np.isnan(g_vs_z)] = 0 - - -# somehow missing a factor of 100 somewhere (with Beer-Lambert)... -# g_vs_z = 100*g_vs_z - -# result: entries are v, efn, and efp. Store in arrays instead. - -# equilibrium solution: -j0, result0 = sesame.IVcurve(sys, [0]) - -voltages = np.concatenate((np.linspace(0,0.55, 8, endpoint=False), - np.linspace(0.55, 0.7, 12, endpoint=False), - np.linspace(0.7, 0.8, 12))) -sys.generation(g_vs_z) -j, result = sesame.IVcurve(sys, voltages)#, verbose=False) -j = j * sys.scaling.current - -jsc = j[0]*1e4 - -# jsc = q*np.trapz(eqe*light_source.spectrum()[1], wls) # A/m2 - -zero_crossing = np.where(np.diff(np.sign(j)))[0][0] -j_above = j[zero_crossing] -j_below = j[zero_crossing + 1] - -voc = voltages[zero_crossing] + (voltages[zero_crossing + 1] - voltages[zero_crossing]) * j_above / (j_above - j_below) - -vmpp = voltages[np.nanargmax(j*voltages)] -jmpp = j[np.nanargmax(j*voltages)]*1e4 - -FF = vmpp*jmpp/(jsc*voc) - -plt.figure() - -plt.plot(voltages, 1e4*j, '-k', label='Current') -plt.plot(voltages, 1e4*j*voltages, '--r', label='Power') - -plt.xlabel('Voltage [V]') -plt.ylabel(r'Current (A/m$^2$) / Power (W)') -plt.text(0.02, 0.9*jsc, r'Jsc = {:.2f} mA/cm$^2$'.format(jsc/10)) -plt.text(0.02, 0.8*jsc, 'Voc = {:.3f} V'.format(voc)) -plt.text(0.02, 0.7*jsc, 'FF = {:.3f} %'.format(FF*100)) -plt.text(0.02, 0.6*jsc, 'Pmax = {:.3f} W'.format(jmpp*vmpp)) -plt.text(0.02, 0.5*jsc, 'Eff = {:.3f} %'.format(jmpp*vmpp/10)) -plt.text(0.02, 0.4*jsc, r'Jmpp = {:.3f} mA/cm$^2$'.format(jmpp/10)) -plt.text(0.02, 0.3*jsc, 'Vmpp = {:.3f} V'.format(vmpp)) - -plt.ylim(-100, 1.1e4*j[0]) -plt.grid() # add grid -plt.legend() -plt.show() - -# print("Jsc", jsc) -# print("Voc", voc) -# print("eff", np.max(1e4*j*voltages)/1000) -# print("FF", FF) -# print("Vmpp", vmpp) -# print("Jmpp", jmpp) - -line = ((0, 0), (L, 0)) - -# sys, result = sesame.load_sim(f'1dhomo_V_0.gzip') -analyzer = sesame.Analyzer(sys, {key: result[key][0, :] for key in result.keys()}) -x, s = analyzer.line(sys, *line) - -electron_density_eq = analyzer.electron_density(location=line)*sys.scaling.density -hole_density_eq = analyzer.hole_density(location=line)*sys.scaling.density - -result_maxPP = {key: result[key][np.nanargmax(j*voltages), :] for key in result.keys()} -# sys, result = sesame.load_sim(f'1dhomo_V_0.gzip') -analyzer = sesame.Analyzer(sys, result_maxPP) -x, s = analyzer.line(sys, *line) - -electron_density = analyzer.electron_density(location=line) -hole_density = analyzer.hole_density(location=line) - -srh = analyzer.bulk_srh_rr(location=line) -auger = analyzer.auger_rr(location=line) -rad = analyzer.radiative_rr(location=line) - -fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True) -ax1.semilogy(sys.xpts[s]*1e4, electron_density*sys.scaling.density, label='e-') -ax1.semilogy(sys.xpts[s]*1e4, hole_density*sys.scaling.density, label='h+') -ax1.set_xlim(0, 10) -ax1.set_xlabel('x (um)') -ax1.set_ylabel('density (cm-3)') -ax1.legend() - -ax2.semilogy(sys.xpts[s]*1e4, electron_density*sys.scaling.density, label='e-') -ax2.semilogy(sys.xpts[s]*1e4, hole_density*sys.scaling.density, label='h+') -# ax2.plot(x, sys.rho*sys.scaling.density) -ax2.set_xlim(d_bulk*1e6 - 8, d_bulk*1e6) -ax2.set_xlabel('x (um)') -plt.show() - - -fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True) -ax1.semilogy(sys.xpts[s]*1e4, electron_density*sys.scaling.density, label='e-') -ax1.semilogy(sys.xpts[s]*1e4, hole_density*sys.scaling.density, label='h+') -ax1.set_xlim(0, 10) -ax1.set_xlabel('x (um)') -ax1.set_ylabel('density (cm-3)') -ax1.legend() - -ax2.semilogy(sys.xpts[s]*1e4, electron_density*sys.scaling.density, label='e-') -ax2.semilogy(sys.xpts[s]*1e4, hole_density*sys.scaling.density, label='h+') -# ax2.plot(x, sys.rho*sys.scaling.density) -ax2.set_xlim(d_bulk*1e6 - 8, d_bulk*1e6) -ax2.set_xlabel('x (um)') -plt.show() - -fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(5, 3)) -ax1.semilogy(sys.xpts[s]*1e4, srh, label='SRH') -ax1.semilogy(sys.xpts[s]*1e4, auger, label='Auger') -ax1.semilogy(sys.xpts[s]*1e4, rad, label='Radiative') -ax1.set_xlim(0, 2) -ax1.set_xlabel('x (um)') -ax1.legend() - -ax2.semilogy(sys.xpts[s]*1e4, srh, label='SRH') -ax2.semilogy(sys.xpts[s]*1e4, auger, label='Auger') -ax2.semilogy(sys.xpts[s]*1e4, rad, label='Radiative') -# ax2.plot(x, sys.rho*sys.scaling.density) -ax2.set_xlim(d_bulk*1e6 - 2, d_bulk*1e6) -ax2.set_xlabel('x (um)') -plt.tight_layout() -plt.show() - -# analyzer.full_current()*sys.scaling.current * sys.scaling.length - -# jsc = analyzer.full_current()*sys.scaling.current - -# bulk SRH: -bulk_SRH = analyzer.integrated_bulk_srh_recombination()*sys.scaling.generation * sys.scaling.length -bulk_aug = analyzer.integrated_auger_recombination()*sys.scaling.generation * sys.scaling.length -bulk_rad = analyzer.integrated_radiative_recombination()*sys.scaling.generation * sys.scaling.length - -# current_loss = q*(bulk_SRH + bulk_aug + bulk_rad)/jsc - -# integrated current: - -jsc_max = 1e9*q*np.trapz(A_intgr*base_light.spectrum(wls*1e9)[1], wls) - -# current_loss_2 = 1 - jsc*1e4/jsc_max -# -# print(current_loss, current_loss_2) - -# account for all the losses - -auger_loss = 1e4*q*bulk_aug -srh_loss = 1e4*q*bulk_SRH -rad_loss = 1e4*q*bulk_rad - -# current available above 1200 nm: - -total_J = 1e9*q*np.trapz(base_light.spectrum(wls*1e9)[1], wls) - -R_loss = 1e9*q*np.trapz(tmm_result['R']*base_light.spectrum(wls*1e9)[1], wls) - -escape_loss = tmm_result['T'] -escape_loss = 1e9*q*np.trapz(escape_loss*base_light.spectrum(wls*1e9)[1], wls) - -parasitic_abs = np.sum(tmm_result['A_per_layer'], 1) - A - -parasitic_loss = 1e9*q*np.trapz(parasitic_abs*base_light.spectrum(wls*1e9)[1], wls) -# units of SRV: cm/s -# units of carrier density: cm-3 -# SRV*carrier density = cm-2/s -# q*SRV*carrier density = q*cm-2/s = A cm-2 - -shading_loss = 1e9*q*np.trapz(shading*base_light.spectrum(wls*1e9)[1], wls) - -# think this should be excess carrier density -n_front_eq = electron_density_eq[0] -p_rear_eq = hole_density_eq[-1] - -n_rear_eq = electron_density_eq[-1] -p_front_eq = hole_density_eq[0] - -e_sr_left = Sn_left*1e4*q*(electron_density[0]*sys.scaling.density - n_front_eq) -h_sr_left = Sp_left*1e4*q*(hole_density[0]*sys.scaling.density - p_front_eq) -e_sr_right = Sn_right*1e4*q*(electron_density[-1]*sys.scaling.density - n_rear_eq) -h_sr_right = Sp_right*1e4*q*(hole_density[-1]*sys.scaling.density - p_rear_eq) - -# only minority carrier recombination matters -sr_left = e_sr_left if hole_density[0]*sys.scaling.density > electron_density[0]*sys.scaling.density else h_sr_left -sr_right = e_sr_right if hole_density[-1]*sys.scaling.density > electron_density[-1]*sys.scaling.density else h_sr_right - -print(shading_loss + R_loss + parasitic_loss + escape_loss + jmpp + auger_loss + srh_loss + rad_loss + sr_left + sr_right) -print(total_J) - -J_list = [jmpp, shading_loss, R_loss, parasitic_loss, escape_loss, auger_loss, srh_loss, rad_loss, sr_left, sr_right] -fig, ax = plt.subplots(figsize=(4, 6)) - -labels = [r'$J_{MPP}$', 'Shading', r'R$_{front}$', r'A$_{parasitic}$', 'Escape', 'Auger', 'SRH', 'Radiative', 'Front surface', 'Rear surface'] - -for j1, item in enumerate(J_list): - - ax.bar(['J'], item, bottom=np.sum(J_list[:j1]), label=labels[j1]) - -ax.set_ylim(300, 1.01*total_J) -ax.axhline(total_J, color='k', linestyle='--', label='Available current') -plt.legend() -plt.show() - -print(total_J - np.sum(J_list)) - -# LONGi: -# Jsc: 41.16 -# Voc: 0.751 -# FF: 86.5 - -# 33.91 -# 0.777 -# 22.766 \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 7d646f3d..cd9dca70 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ [build-system] build-backend = 'mesonpy' requires = [ - 'meson-python==0.13.0rc0', + 'meson-python', 'numpy', 'cython' ] @@ -39,7 +39,7 @@ dependencies = [ "pyyaml", "yabox", "joblib", - "sesame@git+https://github.com/phoebe-p/sesame.git" + "solsesame", ] dynamic = ["version"] @@ -64,12 +64,12 @@ doc = ["Sphinx", "recommonmark"] dev = ["pre-commit"] -[tool.devpy] +[tool.spin] package = 'solcore' -[tool.devpy.commands] -"Build" = ["devpy.cmds.meson.build", "devpy.cmds.meson.test"] -"Extensions" = ['.devpy/cmds.py:codecov', '.devpy/cmds.py:install_dependencies'] +[tool.spin.commands] +"Build" = ["spin.cmds.meson.build", "spin.cmds.meson.test"] +"Extensions" = ['.spin/cmds.py:codecov', '.spin/cmds.py:install_dependencies'] [tool.pytest.ini_options] addopts = "--cov=solcore --cov-report=html:htmlcov -p no:warnings -n \"auto\" -v" diff --git a/solcore/analytic_solar_cells/depletion_approximation.py b/solcore/analytic_solar_cells/depletion_approximation.py index 05ce76e1..37bbef4f 100755 --- a/solcore/analytic_solar_cells/depletion_approximation.py +++ b/solcore/analytic_solar_cells/depletion_approximation.py @@ -364,9 +364,12 @@ def get_j_dark(x, w, L, s, d, V, minority, T): :return: J_top_dark """ - # We calculate some fractions harg = (x - w) / L + + # only do cosh and sinh if harg is < 200 (avoid inf/nan errors) + harg[harg > 200] = 200 + sinh_harg = np.sinh(harg) cosh_harg = np.cosh(harg) lsod = (L * s) / d diff --git a/solcore/material_data/mobility.py b/solcore/material_data/mobility.py index 61675054..59ffd61f 100755 --- a/solcore/material_data/mobility.py +++ b/solcore/material_data/mobility.py @@ -48,7 +48,7 @@ def calculate_mobility(material, holes, N, x=0.0, y=0.0, T=300): if holes: i = 2 if material not in data.keys(): - print("Error: Material {0} not in the database for the mobility. Reverting to GaAs.".format(material)) + print("Warning: material {0} not in the database for the mobility. Reverting to GaAs.".format(material)) d = data['GaAs'][i] elif data[material][0] == 2: d = data[material][i] diff --git a/solcore/material_system/material_system.py b/solcore/material_system/material_system.py index b3ab1817..1f3799ac 100755 --- a/solcore/material_system/material_system.py +++ b/solcore/material_system/material_system.py @@ -353,30 +353,6 @@ def __getattr__(self, attrname): # only used for unknown attributes. return 1.0 / self.ni ** 2 * 2 * pi / (h ** 3 * c ** 2) * quad(inter, 0, upper)[0] if attrname == "permittivity": return self.relative_permittivity * vacuum_permittivity - if attrname == "electron_diffusion_length": - try: - mobility, lifetime = self.electron_mobility, self.electron_minority_lifetime - return np.sqrt(kb * self.T * mobility * lifetime / q) - except: - ValueError("Cannot calculate electron diffusion length. Need electron mobility and lifetime") - if attrname == "hole_diffusion_length": - try: - mobility, lifetime = self.hole_mobility, self.hole_minority_lifetime - return np.sqrt(kb * self.T * mobility * lifetime / q) - except: - ValueError("Cannot calculate hole diffusion length. Need hole mobility and lifetime") - if attrname == "electron_minority_lifetime": - try: - mobility, diffusion_length = self.electron_mobility, self.electron_diffusion_length - return q*diffusion_length**2/(kb*self.T*mobility) - except: - ValueError("Cannot calculate electron minority lifetime. Need electron mobility and diffusion length") - if attrname == "hole_minority_lifetime": - try: - mobility, diffusion_length = self.hole_mobility, self.hole_diffusion_length - return q * diffusion_length ** 2 / (kb * self.T * mobility) - except: - ValueError("Cannot calculate hole minority lifetime. Need hole mobility and diffusion length") kwargs = {element: getattr(self, element) for element in self.composition} kwargs["T"] = self.T diff --git a/solcore/parameter_system/parameter_system.py b/solcore/parameter_system/parameter_system.py index 73be1685..4ff53a1b 100755 --- a/solcore/parameter_system/parameter_system.py +++ b/solcore/parameter_system/parameter_system.py @@ -72,8 +72,8 @@ def get_parameter(self, material, parameter, verbose=False, **others): """ - warn('Passing alloy fractions to get_parameter in the material name (e.g. In0.2GaAs) is now deprecated; \ - pass the alloy fraction as an argument instead.', DeprecationWarning) + # warn('Passing alloy fractions to get_parameter in the material name (e.g. In0.2GaAs) is now deprecated; \ + # pass the alloy fraction as an argument instead.', DeprecationWarning) relevant_parameters = others diff --git a/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py b/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py index 568bff0d..3667559a 100644 --- a/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py +++ b/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py @@ -21,7 +21,7 @@ reason_to_exclude = None except ImportError: reason_to_exclude = ( - "The built-in Poisson-Drift-Difussion solver could not be found." + "The built-in Fortran Poisson-Drift-Diffusion solver could not be found." ) @@ -295,6 +295,10 @@ def iv_pdd( This is in addition to the attributes added by the equilibrium and the short circuit calculations, as relevant. See the description of those functions. """ + + if np.max(internal_voltages) <= 0: + internal_voltages = np.append(internal_voltages, 0.01) + junction.voltage = internal_voltages R_shunt = min(junction.R_shunt, 1e14) if hasattr(junction, "R_shunt") else 1e14 diff --git a/solcore/sesame_drift_diffusion/__init__.py b/solcore/sesame_drift_diffusion/__init__.py index 61c2cc9b..c4b27a7f 100644 --- a/solcore/sesame_drift_diffusion/__init__.py +++ b/solcore/sesame_drift_diffusion/__init__.py @@ -1 +1 @@ -from .solve_pdd import qe_sesame, iv_sesame, equilibrium \ No newline at end of file +from .solve_pdd import qe_sesame, iv_sesame, equilibrium diff --git a/solcore/sesame_drift_diffusion/process_structure.py b/solcore/sesame_drift_diffusion/process_structure.py new file mode 100644 index 00000000..926e6352 --- /dev/null +++ b/solcore/sesame_drift_diffusion/process_structure.py @@ -0,0 +1,522 @@ +from __future__ import annotations + +from sesame import Builder +import numpy as np +from scipy.optimize import root +from solcore.constants import q, kb +from solcore.state import State +from solcore.structure import Junction +from solcore import material + +def process_structure(junction: Junction, options: State): + """ + Process the layers in the junction and convert it to a format which can be used by Sesame. This includes unit + conversions from Solcore base SI units to the units used by Sesame (eV, cm). This function also sets the mesh, + either from a user-supplied junction.mesh (in m) attribute, or by trying to guess an appropriate mesh depending on the + doping profile and layer thicknesses. It will also extract the doping profile, which can be supplied in + the following formats: + + - a constant value for each layer, set by the Nd or Na attribute of the material + - a doping profile for the whole junction, passed as a function which accepts a depth in the junction in m + and returns doping in m-3 + - a doping profile for each layer, passed as a function which accepts a depth in the layer in m and returns + doping in m-3 + - Combinations of the above + + Note that since we can specify a doping profile for the whole junction, the Sesame solver can support a junction + which is made up of a single layer, unlike the depletion approximation and Fortran-based PDD solvers. + + This function will add or update the following attributes in the Junction object: + + - mesh: the mesh used by the Sesame solver, in m + - mesh_cm: the mesh used by the Sesame solver, in cm + - sesame_sys: a Sesame Builder object defining the junction in Sesame format, used to solve the IV and QE + + :param junction: a Junction object + :param options: a State object containing options for the solver + + """ + + materials = [] # list of dictionaries for Sesame + layer_widths = [] # layer widths in cm + doping_profile_functions = ( + [] + ) # list of functions which take an argument of depth in junction in m and + # return doping in cm-3 + offset = 0 # offset relative to front of junction of the current layer in cm + + for layer in junction: + layer_widths.append(layer.width * 1e2) # m to cm + + # this will always be set, even if the layer has a doping profile + # constant_doping = layer.material.Nd / 1e6 if layer.material.Nd > layer.material.Na else -layer.material.Na /1e6 + + if hasattr(layer, "doping_profile"): + doping_profile_functions.append( + lambda i, x: junction[i].doping_profile(x) / 1e6 + ) + # function which takes an argument of depth in junction in m and returns doping in cm-3 + + else: + doping_profile_functions.append( + lambda i, x: junction[i].material.Nd / 1e6 + if junction[i].material.Nd > junction[i].material.Na + else -junction[i].material.Na / 1e6 + ) + # to be consistent with functions passed by user, this should be a function which takes an argument in m + # and returns doping in cm-3 + + offset += layer.width * 1e2 + # Note: in Sesame, all of these can be functions of position as well as constants - but this is untested! + + new_mat = get_material_parameters(layer.material) + materials.append(new_mat) + + # see if we can determine the junction depth: + + edges = np.insert( + np.cumsum(layer_widths), 0, 0 + ) # locations of interfaces between layers + edges[-1] = ( + edges[-1] + 1e-10 + ) # otherwise final point will not be assigned any values + + if hasattr(junction, "doping_profile"): + junction_depth = ( + 100 * root(junction.doping_profile, np.sum(layer_widths) / (2 * 100)).x + ) # in cm + # this is where doping crosses zero - not necessarily the junction depth, but not sure how best to determine + # this. Maybe want to find slope of doping profile and identify locations where it is steepest? + + else: + # identify where the doping profile changes sign for constant doping layers or individual doping profiles + # per layer + doping_sign = [ + np.sign(doping_profile_functions[i1](i1, edges[i1] / 100)) + for i1 in range(len(edges) - 1) + ] + junction_depth = edges[np.where(np.diff(doping_sign))[0] + 1] + + if not hasattr(junction, "mesh"): + # calculate appropriate mesh (hopefully) if not supplied as input in junction.mesh + make_mesh(junction, layer_widths, options, junction_depth) + user_mesh = False + + else: + # user-provided mesh + junction.mesh_cm = junction.mesh * 1e2 # m to cm + user_mesh = True + + if hasattr(junction, "doping_profile"): + doping_profile_x = junction.doping_profile(junction.mesh) / 1e6 # in cm-3 + + else: + doping_profile_x = np.zeros(len(junction.mesh_cm)) + for i1, doping_profile_fn in enumerate(doping_profile_functions): + doping_profile_x[ + np.all( + (junction.mesh_cm >= edges[i1], junction.mesh_cm < edges[i1 + 1]), + axis=0, + ) + ] = doping_profile_fn( + i1, + junction.mesh[ + np.all( + ( + junction.mesh_cm >= edges[i1], + junction.mesh_cm < edges[i1 + 1], + ), + axis=0, + ) + ] + - edges[i1] / 100, + ) + + # Make Sesame Builder object for simulations + junction.sesame_sys = Builder(junction.mesh_cm, T=options.T) # Sesame system + + junction.sesame_sys.rho = doping_profile_x / junction.sesame_sys.scaling.density + + if ( + not user_mesh + ): # do not update mesh after finding doping profile if a specific mesh was supplied by the user + update_mesh(junction, layer_widths, options) + + if hasattr(junction, "doping_profile"): + doping_profile_x = junction.doping_profile(junction.mesh) / 1e6 # in cm-3 + + else: + doping_profile_x = np.zeros(len(junction.mesh_cm)) + for i1, doping_profile_fn in enumerate(doping_profile_functions): + doping_profile_x[ + np.all( + ( + junction.mesh_cm >= edges[i1], + junction.mesh_cm < edges[i1 + 1], + ), + axis=0, + ) + ] = doping_profile_fn( + i1, + junction.mesh[ + np.all( + ( + junction.mesh_cm >= edges[i1], + junction.mesh_cm < edges[i1 + 1], + ), + axis=0, + ) + ] + - edges[i1] / 100, + ) + + junction.sesame_sys = Builder(junction.mesh_cm, T=options.T) # Sesame system + + # set doping profile in Sesame Builder objects + junction.sesame_sys.rho = doping_profile_x / junction.sesame_sys.scaling.density + + for i1 in range(len(junction)): + junction.sesame_sys.add_material( + materials[i1], lambda x: np.all((x >= edges[i1], x < edges[i1 + 1]), axis=0) + ) + + junction.sesame_sys.contact_type("Ohmic", "Ohmic") + # should also be able to choose other options, since Sesame can handle them + + # get surface recombination velocities + junction.sesame_sys.contact_S(*get_srv(junction)) + +def get_material_parameters(mat: material): + """ + Extract approriate material parameters for Sesame from a Solcore material object. + + :param mat: a Solcore material + """ + try: + bulk_recombination_energy = mat.bulk_recombination_energy / q + # this is an option which can be used by Sesame, but is not currently used by Solcore + # anywhere else, or defined in the database for any materials + except ValueError: + bulk_recombination_energy = 0 # default value + + try: + auger_electron = ( + mat.electron_auger_recombination * 1e12 + ) # m6/s to cm6/s + + except ValueError: + auger_electron = 0 + + try: + auger_hole = mat.hole_auger_recombination * 1e12 # m6/s to cm6/s + + except ValueError: + auger_hole = 0 + + # mobility must be provided, but if the user has provided a diffusion length, we can calculate + # the minority carrier lifetime from the mobility and diffusion length + try: + electron_minority_lifetime = mat.electron_minority_lifetime + except ValueError: + electron_minority_lifetime = carrier_constants("electron_minority_lifetime", mat) + + try: + hole_minority_lifetime = mat.hole_minority_lifetime + except ValueError: + hole_minority_lifetime = carrier_constants("hole_minority_lifetime", mat) + + new_mat = { + "Nc": mat.Nc + * 1e-6, # effective density of states at CB edge (cm-3) + "Nv": mat.Nv + * 1e-6, # effective density of states at VB edge (cm-3) + "Eg": mat.band_gap / q, # material bandgap (eV) + "affinity": mat.electron_affinity / q, # electron affinity (eV) + "epsilon": mat.relative_permittivity, + # relative permittivity (dimensionless) + "mu_e": mat.electron_mobility + * 1e4, # electron mobility (cm2/(V s)) + "mu_h": mat.hole_mobility * 1e4, # hole mobility (cm2/(V s)) + "tau_e": electron_minority_lifetime, # electron bulk lifetime (s) + "tau_h": hole_minority_lifetime, # hole bulk lifetime (s) + "Et": bulk_recombination_energy, # energy level of bulk recombination centres (eV) + "B": mat.radiative_recombination * 1e6, + # radiative recombination constant (cm3/s) + "Cn": auger_electron, # Auger recombination constant for electrons (cm6/s) + "Cp": auger_hole, # Auger recombination constant for holes (cm6/s) + } + + return new_mat + +def carrier_constants(quantity, mat): + if quantity == "electron_diffusion_length": + mobility, lifetime = mat.electron_mobility, mat.electron_minority_lifetime + return np.sqrt(kb * mat.T * mobility * lifetime / q) + + if quantity == "hole_diffusion_length": + mobility, lifetime = mat.hole_mobility, mat.hole_minority_lifetime + return np.sqrt(kb * mat.T * mobility * lifetime / q) + + if quantity == "electron_minority_lifetime": + mobility, diffusion_length = mat.electron_mobility, mat.electron_diffusion_length + return q * diffusion_length ** 2 / (kb * mat.T * mobility) + + if quantity == "hole_minority_lifetime": + mobility, diffusion_length = mat.hole_mobility, mat.hole_diffusion_length + return q * diffusion_length ** 2 / (kb * mat.T * mobility) + + +def make_mesh( + junction: Junction, layer_width: float, options: State, junction_depth: list[float] +): + """ + Make a mesh for Sesame, in case one was not supplied. This will try to make the mesh dense near + the front surface, at the (guessed) junction depth, and near the rear surface. The user can specify + the minimum and maximum spacing allowed (in m) through options.minimum_spacing or options.maximum_spacing. + Otherwise these will be set so that there is a maximum of 10,000 and a minimum of 500 points in the + mesh + + :param junction: a Junction object + :param layer_width: a list of layer widths in m + :param options: a State object containing options for the solver + :param junction_depth: a list of junction depths in cm + """ + + # TODO: general improvement of the meshing; should generate a smooth mesh rather than either coarse or fine + + total_width = np.sum(layer_width) + + if "minimum_spacing" in options.keys(): + minimum_spacing = options.minimum_spacing * 1e2 + + else: + minimum_spacing = total_width / 1e4 # maximum of 10000 points + + if "maximum_spacing" in options.keys(): + maximum_spacing = options.maximum_spacing * 1e2 + + else: + maximum_spacing = total_width / 500 # minimum of 500 points + + if maximum_spacing == minimum_spacing: + mesh = np.arange(0, total_width, maximum_spacing) + + # identify location of junction. doping_profile is a function + + else: + # edges = np.insert(np.cumsum(layer_width), 0, 0) + + front_spacing = np.min([minimum_spacing, 0.5e-7]) + + dense_front_width = 200e-7 + + end_front = np.min([dense_front_width - 1e-10, total_width]) + + front = np.arange(0, end_front, front_spacing) + + # 1 nm in cm = 1e-7 + if end_front != total_width: + if junction_depth[0] < total_width / 2: + end_dense = ( + junction_depth[0] + 0.1 * total_width + if (junction_depth[0] + 0.1 * total_width < total_width) + else total_width + ) + + include_end = True if end_dense == total_width else False + + if end_dense > end_front: + n_points = int((end_dense - end_front) / minimum_spacing) + mesh = np.linspace( + end_front, end_dense, n_points, endpoint=include_end + ) + + mesh = np.concatenate((front, mesh)) + + else: + mesh = front + + if not include_end: + n_points = int((total_width - end_dense) / maximum_spacing) + mesh_rest = np.linspace( + end_dense, total_width, n_points, endpoint=True + ) + + mesh = np.concatenate((mesh, mesh_rest)) + + else: + n_points = int((0.2 * total_width - end_front) / minimum_spacing) + + mesh_dense_front = np.linspace( + end_front, 0.2 * total_width, n_points, endpoint=False + ) + + n_points = int( + (junction_depth[0] - 0.3 * total_width) / maximum_spacing + ) + + mesh_bulk = np.linspace( + 0.2 * total_width, + junction_depth[0] - 0.1 * total_width, + n_points, + endpoint=False, + ) + + end_dense = ( + junction_depth[0] + 0.1 * total_width + if (junction_depth[0] + 0.1 * total_width < total_width) + else total_width + ) + + include_end = True if end_dense == total_width else False + + n_points = int( + (end_dense - (junction_depth[0] - 0.1 * total_width)) + / minimum_spacing + ) + + mesh_dense = np.linspace( + junction_depth[0] - 0.1 * total_width, + end_dense, + n_points, + endpoint=include_end, + ) + + if not include_end: + n_points = int((total_width - end_dense) / maximum_spacing) + + mesh_rest = np.linspace( + end_dense, total_width, n_points, endpoint=True + ) + + mesh = np.unique( + np.concatenate( + (front, mesh_dense_front, mesh_bulk, mesh_dense, mesh_rest) + ) + ) + + else: + mesh = np.unique( + np.concatenate((front, mesh_dense_front, mesh_bulk, mesh_dense)) + ) + + junction.mesh_cm = mesh + junction.mesh = junction.mesh_cm * 1e-2 # cm to m + + +def update_mesh(junction: Junction, layer_width: list[float], options: State): + """ + Update mesh after the doping profile is known. + + :param junction: a Junction object + :param layer_width: a list of layer widths in cm + :param options: a State object containing options for the solver + """ + + total_width = np.sum(layer_width) + + # always do 1 nm spacing near front surface? + + if "minimum_spacing" in options.keys(): + minimum_spacing = options.minimum_spacing * 1e2 + + else: + minimum_spacing = total_width / 1e4 # maximum of 10000 points + + current_mesh = junction.mesh_cm + + diff_rho = np.abs(np.gradient(junction.sesame_sys.rho, current_mesh)) + + doping_change = np.where(diff_rho / np.max(diff_rho) > 1e-4)[0] + + # identify continuous stretches of high doping change, and where there are gaps: + + change_point = np.where(np.diff(doping_change) != 1)[0] + + last_points = doping_change[change_point] + first_points = doping_change[change_point + 1] + + first_points = np.insert(first_points, 0, doping_change[0]) + last_points = np.append(last_points, doping_change[-1]) + + first_x = current_mesh[first_points] + last_x = current_mesh[last_points] + + for i1 in range(len(first_points)): + # want to make mesh denser where there are large changes in doping + # print('denser mesh between', first_x[i1], last_x[i1]) + + mesh_above = current_mesh[current_mesh < first_x[i1]] + mesh_below = current_mesh[current_mesh >= last_x[i1]] + mesh_dense = np.arange( + first_x[i1], last_x[i1] - 1e-12, minimum_spacing + ) # last point not included + + current_mesh = np.unique(np.concatenate((mesh_above, mesh_dense, mesh_below))) + + junction.mesh_cm = current_mesh + junction.mesh = junction.mesh_cm * 1e-2 # cm to m + + +def get_srv(junction: Junction): + """ + Extract approriate surface recombination velocities for Sesame from a Junction object. + The labelling sn and sp in the junction object refers to the doping of the region, NOT the + carrier type. For the Sesame solver, it is also possible to specify the recombination velocities + for the carrier types separately in the defintion of the Junction as sn_e/sn_h, sp_e/sp_h, + being respectively the electron/hole SRVs in the n-type region and the electron/hole SRVs in the + p_type region. Note that the majority carrier surface recombination (sn_e and sp_h) generally + do not affect the solution. + + :param junction: a Junction object + """ + + if junction.sesame_sys.rho[junction.sesame_sys.nx - 1] < 0: + # p-type at back + polarity = "np" + + else: + polarity = "pn" + + if hasattr(junction, "sn_e"): + Sn_e = junction.sn_e * 100 + + elif hasattr(junction, "sn"): + Sn_e = junction.sn * 100 + + else: + Sn_e = 1e4 # same as Fortran PDD default (this one is in cm s-1, not m s-1 + + if hasattr(junction, "sp_e"): + Sp_e = junction.sp_e * 100 + + elif hasattr(junction, "sp"): + Sp_e = junction.sp * 100 + + else: + Sp_e = 1e4 + + if hasattr(junction, "sn_h"): + Sn_h = junction.sn_h * 100 + + elif hasattr(junction, "sn"): + Sn_h = junction.sn * 100 + + else: + Sn_h = 1e4 # is this the same as the Fortran PDD? + + if hasattr(junction, "sp_h"): + Sp_h = junction.sp_h * 100 + + elif hasattr(junction, "sp"): + Sp_h = junction.sp * 100 + + else: + Sp_h = 1e4 + + if polarity == "pn": + return Sp_e, Sp_h, Sn_e, Sn_h + + else: + return Sn_e, Sn_h, Sp_e, Sp_h + diff --git a/solcore/sesame_drift_diffusion/solve_pdd.py b/solcore/sesame_drift_diffusion/solve_pdd.py index b83ce48e..b6add5f7 100644 --- a/solcore/sesame_drift_diffusion/solve_pdd.py +++ b/solcore/sesame_drift_diffusion/solve_pdd.py @@ -1,388 +1,79 @@ +from __future__ import annotations + from sesame import Builder, IVcurve, Analyzer import numpy as np from scipy.optimize import root -from solcore.constants import q +from solcore.constants import q, kb from scipy.interpolate import interp1d from solcore.state import State +from solcore.structure import Junction, Layer +from solcore.sesame_drift_diffusion.process_structure import process_structure import warnings -from joblib import Parallel, delayed - -from solcore.registries import ( - register_iv_solver, - register_equilibrium_solver -) - -try: - import sesame - reason_to_exclude = None -except ImportError: - reason_to_exclude = ( - "Sesame was not installed." - ) - -def process_structure(junction, options): +from solcore.registries import register_iv_solver, register_equilibrium_solver +@register_equilibrium_solver("sesame_PDD") +def equilibrium(junction: Junction, + **kwargs): """ - Process the layers in the junction and convert it to a format which can be used by Sesame. This includes unit - conversions from Solcore base SI units to the units used by Sesame (eV, cm). This function also sets the mesh, - either from a user-supplied junction.mesh (in m) attribute, or by guessing an appropriate mesh depending on the - doping profile and layer thicknesses. It will also extract the doping profile, which can be supplied in - the following formats: - - - a constant value for each layer, set by the Nd or Na attribute of the material - - a doping profile for the whole junction, passed as a function which accepts a depth in the junction in m - and returns doping in m-3 - - a doping profile for each layer, passed as a function which accepts a depth in the layer in m and returns - doping in m-3 - - Combinations of the above - - Note that since we can specify a doping profile for the whole junction, the Sesame solver can support a junction - which is made up of a single layer, unlike the depletion approximation and Fortran-based PDD solvers. - - This function will add or update the following attributes in the Junction object: - - - mesh: the mesh used by the Sesame solver, in m - - mesh_cm: the mesh used by the Sesame solver, in cm - - sesame_sys: a Sesame Builder object defining the junction in Sesame format, used to solve the IV and QE - - :param junction: a Junction object - :param options: a State object containing options for the solver - - """ - - materials = [] # list of dictionaries for Sesame - layer_widths = [] # layer widths in cm - doping_profile_functions = [] # list of functions which take an argument of depth in junction in m and - # return doping in cm-3 - offset = 0 # offset relative to front of junction of the current layer in cm - - for layer in junction: - - try: - bulk_recombination_energy = layer.material.bulk_recombination_energy - - except: - bulk_recombination_energy = 0 # should this be the degau - - layer_widths.append(layer.width*1e2) # m to cm - - # this will always be set, even if the layer has a doping profile - constant_doping = layer.material.Nd / 1e6 if layer.material.Nd > layer.material.Na else -layer.material.Na /1e6 - - if hasattr(layer, "doping_profile"): - doping_profile_functions.append(lambda x: layer.doping_profile(x) / 1e6) - # function which takes an argument of depth in junction in m and returns doping in cm-3 - - else: - doping_profile_functions.append(interp1d([offset/100, offset + layer.width], - [constant_doping, constant_doping])) - # to be consistent with functions passed by user, this should be a function which takes an argument in m - # and returns doping in cm-3 - - try: - auger_electron = layer.material.electron_auger_recombination * 1e12 # m6/s to cm6/s - - except: - auger_electron = 0 - - try: - auger_hole = layer.material.hole_auger_recombination * 1e12 # m6/s to cm6/s - - except: - auger_hole = 0 - - new_mat = { - 'Nc': layer.material.Nc * 1e-6, # effective density of states at CB edge (cm-3) - 'Nv': layer.material.Nv * 1e-6, # effective density of states at VB edge (cm-3) - 'Eg': layer.material.band_gap / q, # material bandgap (eV) - 'affinity': layer.material.electron_affinity / q, # electron affinity (eV) - 'epsilon': layer.material.relative_permittivity, # relative permittivity (dimensionless) - 'mu_e': layer.material.electron_mobility*1e4, # electron mobility (cm2/(V s)) - 'mu_h': layer.material.hole_mobility*1e4, # hole mobility (cm2/(V s)) - 'tau_e': layer.material.electron_minority_lifetime, # electron bulk lifetime (s) - 'tau_h': layer.material.hole_minority_lifetime, # hole bulk lifetime (s) - 'Et': bulk_recombination_energy, # energy level of bulk recombination centres (eV) - 'B': layer.material.radiative_recombination*1e6, # radiative recombination constant (cm3/s) - 'Cn': auger_electron, # Auger recombination constant for electrons (cm6/s) - 'Cp': auger_hole, # Auger recombination constant for holes (cm6/s) - } - - offset += layer.width*1e2 - # Note: in Sesame, all of these can be functions of position as well as constants - but this is untested! - - materials.append(new_mat) - - # see if we can determine the junction depth: - - edges = np.insert(np.cumsum(layer_widths), 0, 0) # locations of interfaces between layers - edges[-1] = edges[-1] + 1e-10 # otherwise final point will not be assigned any values - - if hasattr(junction, "doping_profile"): - junction_depth = 100*root(junction.doping_profile, np.sum(layer_widths) / (2*100)).x # in cm - # this is where doping crosses zero - not necessarily the junction depth, but not sure how best to determine - # this. Maybe want to find slope of doping profile and identify locations where it is steepest? - - else: - # identify where the doping profile changes sign for constant doping layers or individual doping profiles - # per layer - doping_sign = [np.sign(doping_profile_functions[i1](edges[i1]+1e-9)) for i1 in range(len(edges) - 1)] - junction_depth = edges[np.where(np.diff(doping_sign))[0] + 1] - - if not hasattr(junction, "mesh"): - # calculate appropriate mesh (hopefully) if not supplied as input in junction.mesh - make_mesh(junction, layer_widths, options, junction_depth) - user_mesh = False - - else: - # user-provided mesh - junction.mesh_cm = junction.mesh*1e2 # m to cm - user_mesh = True - - if hasattr(junction, "doping_profile"): - doping_profile_x = junction.doping_profile(junction.mesh) / 1e6 # in cm-3 - - else: - doping_profile_x = np.zeros(len(junction.mesh_cm)) - for i1, doping_profile_fn in enumerate(doping_profile_functions): - doping_profile_x[np.all((junction.mesh_cm >= edges[i1], junction.mesh_cm < edges[i1 + 1]), axis=0)] = ( - doping_profile_fn(junction.mesh[np.all((junction.mesh_cm >= edges[i1], junction.mesh_cm < edges[i1 + 1]), axis=0)])) - - # Make Sesame Builder object for simulations - junction.sesame_sys = Builder(junction.mesh_cm, T=options.T) # Sesame system - - junction.sesame_sys.rho = doping_profile_x / junction.sesame_sys.scaling.density - - if not user_mesh: - update_mesh(junction, layer_widths, options) - - if hasattr(junction, "doping_profile"): - doping_profile_x = junction.doping_profile(junction.mesh) / 1e6 # in cm-3 - - else: - doping_profile_x = np.zeros(len(junction.mesh_cm)) - for i1, doping_profile_fn in enumerate(doping_profile_functions): - doping_profile_x[np.all((junction.mesh_cm >= edges[i1], junction.mesh_cm < edges[i1 + 1]), axis=0)] = ( - doping_profile_fn(junction.mesh[ - np.all((junction.mesh_cm >= edges[i1], junction.mesh_cm < edges[i1 + 1]), - axis=0)])) - - junction.sesame_sys = Builder(junction.mesh_cm, T=options.T) # Sesame system - - junction.sesame_sys.rho = doping_profile_x / junction.sesame_sys.scaling.density - - for i1 in range(len(junction)): - junction.sesame_sys.add_material(materials[i1], - lambda x: np.all((x >= edges[i1], x < edges[i1 + 1]), axis=0)) - - - junction.sesame_sys.contact_type('Ohmic', 'Ohmic') - # should also be able to choose other options, since Sesame can handle them - - # get surface recombination velocities - junction.sesame_sys.contact_S(*get_srv(junction)) - print('mesh points:', len(junction.mesh_cm)) - - -def make_mesh(junction, layer_width, options, junction_depth): - # want the mesh to be dense close to the front surface, where the doping is changing rapidly/at the junction, - # and where material changes happen. - # (changing from n-type to p-type). - # dense linear mesh near surface and junction - # less dense linear mesh in rest of bulk - - total_width = np.sum(layer_width) - - # always do 1 nm spacing near front surface? - - if "minimum_spacing" in options.keys(): - minimum_spacing = options.minimum_spacing*1e2 - - else: - minimum_spacing = total_width / 1e4 # maximum of 10000 points - - if "maximum_spacing" in options.keys(): - maximum_spacing = options.maximum_spacing*1e2 - - else: - maximum_spacing = total_width / 500 # minimum of 500 points - - print(minimum_spacing, maximum_spacing) - - if maximum_spacing == minimum_spacing: - mesh = np.arange(0, total_width, maximum_spacing) - - # identify location of junction. doping_profile is a function - - else: - # edges = np.insert(np.cumsum(layer_width), 0, 0) - - front_spacing = np.min([minimum_spacing, 0.5e-7]) - - dense_front_width = 50e-7 - - end_front = np.min([dense_front_width - 1e-10, total_width]) - - front = np.arange(0, end_front, front_spacing) - - # 1 nm in cm = 1e-7 - if end_front != total_width: - - if junction_depth[0] < total_width / 2: - print('front junction') - end_dense = junction_depth[0] + 0.1 * total_width if ( - junction_depth[0] + 0.1 * total_width < total_width) else total_width - - include_end = True if end_dense == total_width else False - - n_points = int((end_dense - end_front) / minimum_spacing) - mesh = np.linspace(end_front, end_dense, n_points, endpoint=include_end) - - mesh = np.concatenate((front, mesh)) - - if not include_end: - n_points = int((total_width - end_dense) / maximum_spacing) - mesh_rest = np.linspace(end_dense, total_width, n_points, endpoint=True) - - mesh = np.concatenate((mesh, mesh_rest)) - - else: - - n_points = int((0.2 * total_width- end_front) / minimum_spacing) - - mesh_dense_front = np.linspace(end_front, 0.2 * total_width, n_points, endpoint=False) - - n_points = int((junction_depth[0] - 0.3 * total_width) / maximum_spacing) - - mesh_bulk = np.linspace(0.2 * total_width, junction_depth[0] - 0.1 * total_width, n_points, endpoint=False) - - end_dense = junction_depth[0] + 0.1 * total_width if ( - - junction_depth[0] + 0.1 * total_width < total_width) else total_width - - include_end = True if end_dense == total_width else False - - n_points = int((end_dense - (junction_depth[0] - 0.1 * total_width)) / minimum_spacing) - - mesh_dense = np.linspace(junction_depth[0] - 0.1 * total_width, end_dense, n_points, endpoint=include_end) - - if not include_end: - n_points = int((total_width - end_dense) / maximum_spacing) - - mesh_rest = np.linspace(end_dense, total_width, n_points, endpoint=True) + Solve at equilibrium (no illumination, no applied voltage) using the Sesame solver. - mesh = np.unique(np.concatenate((front, mesh_dense_front, mesh_bulk, mesh_dense, mesh_rest))) + :param junction: a Junction object + :param options: a State object containing options for the solver + """ - else: - mesh = np.unique(np.concatenate((front, mesh_dense_front, mesh_bulk, mesh_dense))) + # TODO: pass output from 'result' parameter to the user, similar to the Fortran PDD solver, so that + # band structures etc. can be plotted. - junction.mesh_cm = mesh - junction.mesh = junction.mesh_cm*1e-2 # cm to m + options = State(kwargs) # needs to be passed as kwargs to be compatible with Fortran equilibrium solver + # in registries -def update_mesh(junction, layer_width, options): - - total_width = np.sum(layer_width) - - # always do 1 nm spacing near front surface? - - if "minimum_spacing" in options.keys(): - minimum_spacing = options.minimum_spacing*1e2 - - else: - minimum_spacing = total_width / 1e4 # maximum of 10000 points - - current_mesh = junction.mesh_cm - - diff_rho = np.abs(np.gradient(junction.sesame_sys.rho, current_mesh)) - - doping_change = np.where(diff_rho / np.max(diff_rho) > 1e-4)[0] - - # identify continuous stretches of high doping change, and where there are gaps: - - change_point = np.where(np.diff(doping_change) != 1)[0] - - last_points = doping_change[change_point] - first_points = doping_change[change_point + 1] - - first_points = np.insert(first_points, 0, doping_change[0]) - last_points = np.append(last_points, doping_change[-1]) - - first_x = current_mesh[first_points] - last_x = current_mesh[last_points] - - for i1 in range(len(first_points)): - # want to make mesh denser where there are large changes in doping - # print('denser mesh between', first_x[i1], last_x[i1]) - - mesh_above = current_mesh[current_mesh < first_x[i1]] - mesh_below = current_mesh[current_mesh >= last_x[i1]] - mesh_dense = np.arange(first_x[i1], last_x[i1]-1e-12, minimum_spacing) # last point not included - - current_mesh = np.unique(np.concatenate((mesh_above, mesh_dense, mesh_below))) - - junction.mesh_cm = current_mesh - junction.mesh = junction.mesh_cm*1e-2 # cm to m - -def get_srv(junction): - # Define the surface recombination velocities for electrons and holes [cm/s] - if hasattr(junction, "sn_front"): - Sn_left = junction.sn_front * 100 - - elif hasattr(junction, "sn"): - Sn_left = junction.sn * 100 - - else: - Sn_left = 1e4 # is this the same as the Fortran PDD? - - if hasattr(junction, "sp_front"): - Sp_left = junction.sp_front * 100 - - elif hasattr(junction, "sp"): - Sp_left = junction.sp * 100 - - else: - Sp_left = 1e4 - - if hasattr(junction, "sn_rear"): - Sn_right = junction.sn_rear * 100 + if not hasattr(junction, "sesame_sys"): + process_structure(junction, options) - elif hasattr(junction, "sn"): - Sn_right = junction.sn * 100 + j, result = IVcurve(junction.sesame_sys, [0]) - else: - Sn_right = 1e4 # is this the same as the Fortran PDD? + if np.any(np.isnan(j)): + warnings.warn( + "Current calculation did not converge at all voltages", UserWarning + ) - if hasattr(junction, "sp_rear"): - Sp_right = junction.sp_rear * 100 + j = j * junction.sesame_sys.scaling.current * 1e4 # cm-2 -> m-2 - elif hasattr(junction, "sp"): - Sp_right = junction.sp * 100 + junction.sesame_output = result - else: - Sp_right = 1e4 - return Sn_left, Sp_left, Sn_right, Sp_right +@register_iv_solver("sesame_PDD") +def iv_sesame(junction, options): + """ + Solve the dark or light IV using the Sesame solver. This will scan through the voltages + in options.internal_voltages and call sesame.IVcurve. If your calculation is failing to converge, + make sure your calculation includes 0 V and try scanning more voltage points or setting a denser + mesh. Note that the solver will not calculate anything at voltages > the bandgap + 10*kb*T, since + it will fail to converge for the high injection regime. -@register_equilibrium_solver("sesame_PDD", reason_to_exclude=reason_to_exclude) -def equilibrium(): - pass + :param junction: a Junction object + :param options: a State object containing options for the solver + """ -@register_iv_solver("sesame_PDD", reason_to_exclude=reason_to_exclude) -def iv_sesame(junction, options): + # TODO: pass output from 'result' parameter to the user, similar to the Fortran PDD solver, so that + # band structures etc. can be plotted. if not hasattr(junction, "sesame_sys"): process_structure(junction, options) - # gen_wl = profile_func - - if options.light_iv and np.all(junction.sesame_sys.g == 0): - - gen_wl = junction.absorbed(junction.mesh) / 100 # m-1 -> cm-1 + if options.light_iv: # and np.all(junction.sesame_sys.g == 0): + gen_wl = junction.absorbed(junction.mesh) / 100 # m-1 -> cm-1 wls = options.wavelength - gg = options.light_source.spectrum(wls, output_units="photon_flux_per_m")[1][:, None] * gen_wl.T + gg = ( + options.light_source.spectrum(wls, output_units="photon_flux_per_m")[1][ + :, None + ] + * gen_wl.T + ) - g_vs_z = np.trapz(gg, wls, axis=0) / 1e4 # m^2 -> cm^2 + g_vs_z = np.trapz(gg, wls, axis=0) / 1e4 # m^2 -> cm^2 g_vs_z[np.isnan(g_vs_z)] = 0 # can also pass a function to generation - more flexible? @@ -392,19 +83,32 @@ def iv_sesame(junction, options): R_shunt = min(junction.R_shunt, 1e14) if hasattr(junction, "R_shunt") else 1e14 - # voltages need to go from 0 (or close) to highest applied +ve or -ve voltage. Split if necessary. + max_Eg = np.max(junction.sesame_sys.Eg * junction.sesame_sys.scaling.energy) + max_V = ( + max_Eg + 10 * 8.617e-05 * options.T + ) # do not go into the high injection regime, will not get convergence + + # voltages need to go from 0 (or close) to highest applied +ve or -ve voltage, otherwise + # do not get convergence; need to go from V = 0 to high applied voltage so that Sesame can + # use the previous solution as a guess for the next voltage. - if junction.sesame_sys.rho[junction.sesame_sys.nx-1] < 0: + if junction.sesame_sys.rho[junction.sesame_sys.nx - 1] < 0: # this is necessary because Sesame will internally flip the sign for an n-p junction voltages_for_solve = -voltages if np.all(options.voltages >= 0): - warnings.warn('All voltages are positive, but junction has been identified as n-p, so the ' - 'open-circuit voltage (Voc) of the junction will be negative.', UserWarning) + warnings.warn( + "All voltages are positive, but junction has been identified as n-p, so the " + "open-circuit voltage (Voc) of the junction will be negative.", + UserWarning, + ) else: voltages_for_solve = voltages + voltages_for_solve = voltages_for_solve[voltages_for_solve <= max_V] + + warnings.filterwarnings("ignore") # split +ve and -ve voltages if necessary: if np.any(voltages_for_solve < 0): if np.any(voltages_for_solve > 0): @@ -419,31 +123,43 @@ def iv_sesame(junction, options): negative_voltages = negative_voltages[negative_voltages_order] positive_voltages = positive_voltages[positive_voltages_order] - j_positive, result_positive = sesame.IVcurve(junction.sesame_sys, positive_voltages) - j_negative, result_negative = sesame.IVcurve(junction.sesame_sys, negative_voltages) + j_positive, result_positive = IVcurve( + junction.sesame_sys, positive_voltages + ) + j_negative, result_negative = IVcurve( + junction.sesame_sys, negative_voltages + ) j_negative = j_negative[::-1] - result_positive = {key: result_positive[key][positive_voltages_order, :] for key in result_positive.keys()} - result_negative = {key: result_negative[key][negative_voltages_order, :] for key in result_negative.keys()} + result_negative = { + key: result_negative[key][::-1] + for key in result_negative.keys() + } negative_voltages = negative_voltages[::-1] if np.any(voltages_for_solve == 0): - # V = 0 would have been included in both the +ve and -ve voltages, so # exclude it from the negative voltage results when concatenating j = np.concatenate((j_negative[:-1], j_positive)) - result = {key: np.concatenate((result_negative[key][:-1], result_positive[key])) for key in - result_positive.keys()} - final_voltages = np.concatenate((negative_voltages[:-1], positive_voltages)) - + result = { + key: np.concatenate( + (result_negative[key][:-1], result_positive[key]) + ) + for key in result_positive.keys() + } + final_voltages = np.concatenate( + (negative_voltages[:-1], positive_voltages) + ) else: - j = np.concatenate((j_negative, j_positive)) - result = {key: np.concatenate((result_negative[key], result_positive[key])) for key in result_positive.keys()} + result = { + key: np.concatenate((result_negative[key], result_positive[key])) + for key in result_positive.keys() + } final_voltages = np.concatenate((negative_voltages, positive_voltages)) # this results in j and result in order of increasing values for voltages_for_solve @@ -451,21 +167,29 @@ def iv_sesame(junction, options): # negative voltages only voltage_order = np.argsort(voltages_for_solve)[::-1] final_voltages = voltages_for_solve[voltage_order] - j, result = sesame.IVcurve(junction.sesame_sys, final_voltages) + j, result = IVcurve(junction.sesame_sys, final_voltages) else: # positive voltages only voltage_order = np.argsort(voltages_for_solve) final_voltages = voltages_for_solve[voltage_order] - j, result = sesame.IVcurve(junction.sesame_sys, final_voltages) # , verbose=False) + j, result = IVcurve( + junction.sesame_sys, final_voltages + ) # , verbose=False) + + warnings.resetwarnings() + + if np.any(np.isnan(j)): + warnings.warn( + "Current calculation did not converge at all voltages", UserWarning + ) # final_voltages are the voltages corresponding to the entries in j and result, USING # Sesame's sign convention. So if the voltage sign was flipped above, need to flip it back # for Solcore if junction.sesame_sys.rho[junction.sesame_sys.nx - 1] < 0: - print('flipping back') result_voltage = -final_voltages sort_result = np.argsort(result_voltage) j = j[sort_result] @@ -475,24 +199,47 @@ def iv_sesame(junction, options): else: result_voltage = final_voltages - j = j * junction.sesame_sys.scaling.current * 1e4 # cm-2 -> m-2 + j = j * junction.sesame_sys.scaling.current * 1e4 # cm-2 -> m-2 + + shunted_current = j + result_voltage / R_shunt + + non_nans = np.where(~np.isnan(shunted_current))[0] - junction.current = j + result_voltage / R_shunt + if len(non_nans) > 0: + # may be NaNs in j - find values closest to edge which are not NaN: + first_non_nan = shunted_current[non_nans[0]] + last_non_nan = shunted_current[non_nans[-1]] - # current_increasing = j[voltage_increasing] + else: + Exception( + "No solutions found for IV curve. Try increasing the number of voltage points scanned." + ) + + junction.sesame_output = result junction.iv = interp1d( - result_voltage, - junction.current, + result_voltage[non_nans], + shunted_current[non_nans], kind="linear", bounds_error=False, assume_sorted=True, - fill_value=(j[0], j[-1]), + fill_value=(first_non_nan, last_non_nan), ) - junction.voltage = result_voltage + junction.voltage = options.internal_voltages + junction.current = junction.iv(options.internal_voltages) + junction.pdd_output = process_sesame_results(junction.sesame_sys, result) + + +def qe_sesame(junction: Junction, options: State): + """ + Calculate the quantum efficiency of a junction using Sesame. This will scan through the wavelengths + set in options.wavelength. It will scan from long wavelengths to short wavelengths, to improve + the chance of convergence, since carrier generation profiles will be less steep at longer wavelengths. -def qe_sesame(junction, options): + :param junction: a Junction object + :param options: a State object containing options for the solver + """ if not hasattr(junction, "sesame_sys"): process_structure(junction, options) @@ -501,7 +248,7 @@ def qe_sesame(junction, options): # # else: # n_jobs = options.n_jobs if hasattr(options, "n_jobs") else -1 - n_jobs = 1 + # n_jobs = 1 # parallel implementation does not work due to pickling error. Unsure of cause. wls = options.wavelength @@ -509,48 +256,86 @@ def qe_sesame(junction, options): voltages = [0] # profile_func = interp1d(bulk_positions_cm, 1e7 * Si_profile, kind='linear', bounds_error=False, fill_value=0) - profile_func = junction.absorbed # this returns an array of shape (mesh_points, wavelengths) + profile_func = ( + junction.absorbed + ) # this returns an array of shape (mesh_points, wavelengths) A = np.trapz(junction.absorbed(junction.mesh), junction.mesh, axis=0) def make_gfcn_fun(wl_index, flux): def gcfn_fun(x, y): - return flux * profile_func(x / 100)[wl_index] / 100 # convert to cm-1 from m-1 + return ( + flux * profile_func(np.array([x / 100]))[0, wl_index] / 100 + ) # convert to cm-1 from m-1 return gcfn_fun + # def qe_i(i1): + # + # junction.sesame_sys.generation(make_gfcn_fun(i1, flux)) + # + # j, result = IVcurve(junction.sesame_sys, voltages) + # + # eqe = np.abs(j) / (q * flux) + # + # return eqe, result + # + # allres = Parallel(n_jobs=n_jobs)( + # delayed(qe_i)( + # i1, + # ) + # for i1 in range(len(wls)) + # ) + # + # eqe = np.concatenate([item[0] for item in allres]) + # + # efn = np.concatenate([item[1]['efn'] for item in allres]) + # efp = np.concatenate([item[1]['efp'] for item in allres]) + # vres = np.concatenate([item[1]['v'] for item in allres]) + flux = 1e20 - # for i1, wl in enumerate(wls): + eqe = np.zeros_like(wls) + + # do not solve EQE if absorption is ~ 0 + + EQE_threshold = 1e-5 + + wl_solve = np.where(A >= EQE_threshold)[0][::-1] + # go in backwards order through wavelengths - since generation profile tends to be flatter at longer wavelengths, + # this increases the change of convergence, since the solution for the previous wavelength is always used as a guess + # for the next wavelength. Gaving a good guess helps the short wavelength solutions converge - def qe_i(i1): + warnings.filterwarnings("ignore") + # this is to prevent warning from Sesame flooding the output. Not ideal but unsure on best way to solve. + for i1 in wl_solve: junction.sesame_sys.generation(make_gfcn_fun(i1, flux)) - j, result = sesame.IVcurve(junction.sesame_sys, voltages) + if i1 == wl_solve[0]: + guess = None - eqe = np.abs(j) / (q * flux) + else: + guess = {key: result[key][0, :] for key in result.keys()} - return eqe, result + j, result = IVcurve(junction.sesame_sys, voltages, guess=guess) - allres = Parallel(n_jobs=n_jobs)( - delayed(qe_i)( - i1, - ) - for i1 in range(len(wls)) - ) + eqe[i1] = np.abs(j) / (q * flux) - eqe = np.concatenate([item[0] for item in allres]) + if np.any(np.isnan(eqe)): + warnings.warn( + "EQE calculation did not converge at all wavelengths", UserWarning + ) - efn = np.concatenate([item[1]['efn'] for item in allres]) - efp = np.concatenate([item[1]['efp'] for item in allres]) - vres = np.concatenate([item[1]['v'] for item in allres]) + warnings.resetwarnings() eqe = eqe * junction.sesame_sys.scaling.current - iqe = eqe / A + iqe = np.divide(eqe, A, out=np.zeros_like(eqe), where=A > 0) + + # line = ((0, 0), (np.max(junction.mesh_cm), 0)) + # scale_sesame_result(junction.sesame_sys, result, line) # convert dimensionless current to dimension-ful current - # iqe = iqe * junction.sesame_sys.scaling.current junction.iqe = interp1d(wls, iqe) @@ -563,14 +348,73 @@ def qe_i(i1): fill_value=(eqe[0], eqe[-1]), ) - junction.qe = State( - { - "WL": wls, - "IQE": junction.iqe(wls), - "EQE": junction.eqe(wls), - } - ) + junction.qe = State({"WL": wls, "IQE": junction.iqe(wls), "EQE": junction.eqe(wls)}) +def process_sesame_results(sys: Builder, result: dict): + """ + Scale the result of a Sesame calculation to SI units, and calculate other quantities like + the positions of the conduction and valence band and the recombination rates. Produces a + State object with entries: + + - potential: the potential (V) + - n: the electron density (m\ :sup:`-3`) + - p: the hole density (m\ :sup:`-3`) + - Ec: the level of the conduction band (eV) + - Ev the level of the valence band (eV) + - Efe: the electron quasi-Fermi level (eV) + - Efh the hole quasi-Fermi level (eV) + - Rrad: the radiative recombination rate (m\ :sup:`-3` s\ :sup:`-1`) + - Raug: the Auger recombination rate (m\ :sup:`-3` s\ :sup:`-1`) + - Rsrh: the bulk recombination due to Shockley-Read-Hall processes (m\ :sup:`-3` s\ :sup:`-1`) + +Each of these is a 2-dimensional array, with dimensions ``(len(options.internal_voltages), len(mesh))``. + + + :param sys: a Sesame Builder object + :param result: a dictionary containing the results from a Sesame calculation + """ + + line = ((0, 0), (np.max(sys.xpts), 0)) + n_voltages = len(result["v"]) + + potential = result['v'] * sys.scaling.energy + Efe = result['efn'] * sys.scaling.energy + Efh = result['efp'] * sys.scaling.energy + Ec = -(result['v'] + sys.bl) * sys.scaling.energy + Ev = -(result['v'] + sys.bl + sys.Eg) * sys.scaling.energy + + n = np.zeros((n_voltages, sys.nx)) + p = np.zeros((n_voltages, sys.nx)) + + Rrad = np.zeros((n_voltages, sys.nx)) + Raug = np.zeros((n_voltages, sys.nx)) + Rsrh = np.zeros((n_voltages, sys.nx)) + + for i1 in range(n_voltages): + result_loop = {key: result[key][i1, :] for key in result.keys()} + + analyzer = Analyzer(sys, result_loop) + + n[i1] = analyzer.electron_density(location=line) * sys.scaling.density * 1e6 # m-3 + p[i1] = analyzer.hole_density(location=line) * sys.scaling.density * 1e6 # m-3 + + Rsrh[i1] = analyzer.bulk_srh_rr(location=line) * sys.scaling.generation * 1e6 # m-3 + Raug[i1] = analyzer.auger_rr(location=line) * sys.scaling.generation * 1e6 # m-3 + Rrad[i1] = analyzer.radiative_rr(location=line) * sys.scaling.generation * 1e6 # m-3 + + output = State( + potential=potential, + Efe=Efe, + Efh=Efh, + Ec=Ec, + Ev=Ev, + n=n, + p=p, + Rrad=Rrad, + Raug=Raug, + Rsrh=Rsrh, + ) + return output diff --git a/tests/test_depletion_approximation.py b/tests/test_depletion_approximation.py index 995ce29c..f47babea 100644 --- a/tests/test_depletion_approximation.py +++ b/tests/test_depletion_approximation.py @@ -31,12 +31,16 @@ def test_get_j_dark(): w = wn + wp + xi + arg = (x - w) / L + + arg[arg > 200] = 200 + expected = ( (q * d * minor / L) * (np.exp(q * V / (kb * T)) - 1) * ( - (((s * L) / d) * np.cosh((x - w) / L) + np.sinh((x - w) / L)) - / (((s * L) / d) * np.sinh((x - w) / L) + np.cosh((x - w) / L)) + (((s * L) / d) * np.cosh(arg) + np.sinh(arg)) + / (((s * L) / d) * np.sinh(arg) + np.cosh(arg)) ) ) diff --git a/tests/test_examples.py b/tests/test_examples.py index cfd6ab52..ed782158 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -51,6 +51,8 @@ def test_example_scripts(example, examples_directory): skip("No SMARTS solver found.") elif b"RCWASolverError" in process.stderr: skip("No RCWA solver found.") + elif b"FigureCanvasAgg" in process.stderr: + skip("Plotting issue due to non-interactive Agg backend.") elif process.stderr != b"": raise Exception(process.stderr.decode()) @@ -82,5 +84,7 @@ def test_example_notebooks(example, examples_directory): skip("No SMARTS solver found.") elif b"RCWASolverError" in process.stderr: skip("No RCWA solver found.") + elif b"FigureCanvasAgg" in process.stderr: + skip("Plotting issue due to non-interactive Agg backend.") elif process.stderr != b"": raise Exception(process.stderr.decode()) diff --git a/tests/test_optics.py b/tests/test_optics.py index e34a0379..008657f7 100644 --- a/tests/test_optics.py +++ b/tests/test_optics.py @@ -126,9 +126,9 @@ def test_pol_rcwa(): rat_output_A=A_output_p, parallel=True, steps_size=step_size, pol='p') - assert np.mean(A_output_u, 1) == approx(0.5*(A_output_p + A_output_s)) + assert np.mean(A_output_u, 1) == approx(0.5*(A_output_p + A_output_s), rel=1e-4) - assert result_u["absorption"] == approx(0.5*(result_s["absorption"] + result_p["absorption"])) + assert result_u["absorption"] == approx(0.5*(result_s["absorption"] + result_p["absorption"]), rel=1e-4) @mark.skipif(skip_s4_test(), reason="Only works if S4 installed") @@ -176,7 +176,7 @@ def test_arbitrary_pol_rcwa(): assert A_output == approx(A_output_s) - assert result["absorption"] == approx(result_s["absorption"]) + assert result["absorption"] == approx(result_s["absorption"], rel=1e-4) @mark.skipif(skip_s4_test(), reason="Only works if S4 installed") diff --git a/tests/test_sesame_pdd.py b/tests/test_sesame_pdd.py index fe94cd1d..bc6f8a49 100644 --- a/tests/test_sesame_pdd.py +++ b/tests/test_sesame_pdd.py @@ -1,11 +1,11 @@ import numpy as np -from pytest import approx +from pytest import approx, raises def test_constant_doping(): # test sesame PDD process_structure when constant doping levels are passed from solcore import material, si from solcore.structure import Junction, Layer - from solcore.sesame_drift_diffusion.solve_pdd import process_structure + from solcore.sesame_drift_diffusion.process_structure import process_structure from solcore.state import State Si_n = material('Si')(Nd=1e24, electron_minority_lifetime=1e-6, hole_minority_lifetime=1e-7) @@ -27,28 +27,108 @@ def test_constant_doping(): assert sesame_obj.rho[0] == 1e24*1e-6/sesame_obj.scaling.density assert sesame_obj.rho[-1] == -2e24*1e-6/sesame_obj.scaling.density - pass -def test_layer_doping_profile(): - # test process_structure when a doping profile is passed for one of the layers - pass +def test_doping_profile(): + # test process_structure when a doping profile is passed for one of the layers, + # or to the whole junction + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.sesame_drift_diffusion.process_structure import process_structure + from solcore.state import State + from scipy.interpolate import interp1d + + Si_n = material('Si')(Nd=1e25, electron_minority_lifetime=1e-6, hole_minority_lifetime=1e-7) + Si_i = material('Si')(electron_minority_lifetime=2e-6, hole_minority_lifetime=2e-7) + Si_p = material('Si')(Na=2e24, electron_minority_lifetime=1.5e-6, hole_minority_lifetime=1.5e-7) + + doping_profile = np.linspace(Si_n.Nd, -Si_p.Na, 2000) + x_pos = np.linspace(0, 2000, 2000) * 1e-9 + doping_profile_func = interp1d(x_pos, doping_profile, fill_value=(Si_n.Nd, -Si_p.Na), bounds_error=False) + + junction = Junction([Layer(si('200nm'), Si_n), + Layer(si('2000nm'), Si_i, doping_profile=doping_profile_func), + Layer(si('2000nm'), Si_p)]) + + options = State() + options.T = 300 + + process_structure(junction, options) -def test_junction_doping_profile(): - # test process_structure when a doping profile is passed for the whole junction - pass + rho_1 = interp1d(junction.mesh, junction.sesame_sys.rho) + + x_pos_2 = np.linspace(200, 2200, 2000) * 1e-9 + doping_profile_func = interp1d(x_pos_2, doping_profile, fill_value=(Si_n.Nd, -Si_p.Na), bounds_error=False) + + junction = Junction([Layer(si('200nm'), Si_n), + Layer(si('2000nm'), Si_i), + Layer(si('2000nm'), Si_p)], doping_profile=doping_profile_func) + + process_structure(junction, options) + + rho_2 = interp1d(junction.mesh, junction.sesame_sys.rho) + + assert np.allclose(rho_1(junction.mesh), rho_2(junction.mesh)) def test_parameter_extraction(): # test all material parameters are correctly extracted by process_structure - pass + from solcore.sesame_drift_diffusion.process_structure import get_material_parameters + from solcore import material + from solcore.constants import q, kb + + test_mat = material('Si')(Nd=1e24, Na=2e24) + + with raises(ValueError) as excinfo: + get_material_parameters(test_mat) + assert excinfo.match("nor in calculable parameters") + + test_mat.hole_diffusion_length, test_mat.electron_diffusion_length = 100e-9, 200e-9 + + expected_tau_e = (q * test_mat.electron_diffusion_length ** 2 / + (kb * test_mat.T * test_mat.electron_mobility)) + expected_tau_h = (q * test_mat.hole_diffusion_length ** 2 / + (kb * test_mat.T * test_mat.hole_mobility)) + + result = get_material_parameters(test_mat) + + expected_dict = { + "Nc": test_mat.Nc * 1e-6, + "Nv": test_mat.Nv * 1e-6, + "Eg": test_mat.band_gap / q, + "affinity": test_mat.electron_affinity / q, + "epsilon": test_mat.relative_permittivity, + "mu_e": test_mat.electron_mobility * 1e4, + "mu_h": test_mat.hole_mobility * 1e4, + "tau_e": expected_tau_e, + "tau_h": expected_tau_h, # hole bulk lifetime (s) + "Et": 0, # energy level of bulk recombination centres (eV) + "B": test_mat.radiative_recombination * 1e6, + # radiative recombination constant (cm3/s) + "Cn": test_mat.electron_auger_recombination * 1e12, + # Auger recombination constant for electrons (cm6/s) + "Cp": test_mat.hole_auger_recombination * 1e12, + # Auger recombination constant for holes (cm6/s) + } -def test_dark_IV(): - pass + assert result == expected_dict -def test_light_IV(): - pass + test_mat.electron_minority_lifetime = 1e-6 + test_mat.hole_minority_lifetime = 1e-5 + + expected_dict["tau_e"] = 1e-6 + expected_dict["tau_h"] = 1e-5 + + result = get_material_parameters(test_mat) + + assert result == expected_dict + + test_mat.bulk_recombination_energy = 0.5 + + expected_dict["Et"] = test_mat.bulk_recombination_energy / q + + result = get_material_parameters(test_mat) + + assert result == expected_dict -def test_multijunction_IV(): - pass def test_compare_DA_np(): @@ -56,7 +136,8 @@ def test_compare_DA_np(): from solcore.structure import Junction, Layer from solcore.solar_cell_solver import solar_cell_solver, SolarCell from solcore.analytic_solar_cells import iv_depletion - from solcore.sesame_drift_diffusion import iv_sesame + from solcore.sesame_drift_diffusion.solve_pdd import iv_sesame + from solcore.sesame_drift_diffusion.process_structure import carrier_constants from solcore.state import State from solcore.light_source import LightSource @@ -65,6 +146,9 @@ def test_compare_DA_np(): GaAs_p = material('GaAs')(T=300, Na=1e24, hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_p.electron_diffusion_length = carrier_constants("electron_diffusion_length", GaAs_p) + GaAs_n.hole_diffusion_length = carrier_constants("hole_diffusion_length", GaAs_n) + options = State() options.wavelength = np.linspace(300, 950, 100)*1e-9 options.voltages = np.linspace(-1.3, 0.5, 30) @@ -77,18 +161,18 @@ def test_compare_DA_np(): mesh = np.linspace(0, 2200, 500)*1e-9 - pn_junction = Junction([Layer(si('200nm'), GaAs_n, role='emitter'), Layer(si('2000nm'), GaAs_p, role='base')], - kind='DA', R_shunt=0.1, mesh=mesh) + np_junction = Junction([Layer(si('200nm'), GaAs_n, role='emitter'), Layer(si('2000nm'), GaAs_p, role='base')], + kind='DA', R_shunt=0.1, mesh=mesh, sn=10, sp=10) - solar_cell_solver(SolarCell([pn_junction]), 'optics', options) + solar_cell_solver(SolarCell([np_junction]), 'optics', options) - iv_depletion(pn_junction, options) - depletion_current = pn_junction.current - depletion_current_interp = pn_junction.iv(options.voltages) + iv_depletion(np_junction, options) + depletion_current = np_junction.current + depletion_current_interp = np_junction.iv(options.voltages) - iv_sesame(pn_junction, options) - sesame_current = pn_junction.current - sesame_current_interp = pn_junction.iv(options.voltages) + iv_sesame(np_junction, options) + sesame_current = np_junction.current + sesame_current_interp = np_junction.iv(options.voltages) assert depletion_current[-1] == approx(sesame_current[-1], rel=0.05) assert np.sign(depletion_current[0]) == np.sign(sesame_current[0]) @@ -105,15 +189,19 @@ def test_compare_DA_pn(): from solcore.structure import Junction, Layer from solcore.solar_cell_solver import solar_cell_solver, SolarCell from solcore.analytic_solar_cells import iv_depletion - from solcore.sesame_drift_diffusion import iv_sesame + from solcore.sesame_drift_diffusion.solve_pdd import iv_sesame + from solcore.sesame_drift_diffusion.process_structure import carrier_constants from solcore.state import State from solcore.light_source import LightSource - GaAs_n = material('GaAs')(T=300, Na=1e24, + GaAs_p = material('GaAs')(T=300, Na=1e24, hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) - GaAs_p = material('GaAs')(T=300, Nd=1e24, + GaAs_n = material('GaAs')(T=300, Nd=1e24, hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_p.electron_diffusion_length = carrier_constants("electron_diffusion_length", GaAs_p) + GaAs_n.hole_diffusion_length = carrier_constants("hole_diffusion_length", GaAs_n) + options = State() options.wavelength = np.linspace(300, 950, 100)*1e-9 options.voltages = np.linspace(-0.5, 1.3, 30) @@ -126,7 +214,7 @@ def test_compare_DA_pn(): mesh = np.linspace(0, 2200, 500)*1e-9 - pn_junction = Junction([Layer(si('200nm'), GaAs_n, role='emitter'), Layer(si('2000nm'), GaAs_p, role='base')], + pn_junction = Junction([Layer(si('200nm'), GaAs_p, role='emitter'), Layer(si('2000nm'), GaAs_n, role='base')], kind='DA', R_shunt=0.1, mesh=mesh) solar_cell_solver(SolarCell([pn_junction]), 'optics', options) @@ -146,3 +234,281 @@ def test_compare_DA_pn(): assert np.sign(depletion_current_interp[-1]) == np.sign(sesame_current_interp[-1]) assert np.all(sesame_current == sesame_current_interp) + + +def test_qe(): + + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.solar_cell_solver import solar_cell_solver, SolarCell + from solcore.analytic_solar_cells import qe_depletion + from solcore.sesame_drift_diffusion.solve_pdd import qe_sesame + from solcore.sesame_drift_diffusion.process_structure import carrier_constants + from solcore.state import State + + GaAs_p = material('GaAs')(T=300, Na=1e24, + hole_minority_lifetime=1e-6, + electron_minority_lifetime=1e-6) + GaAs_n = material('GaAs')(T=300, Nd=1e24, + hole_minority_lifetime=1e-6, + electron_minority_lifetime=1e-6) + + GaAs_p.electron_diffusion_length = carrier_constants("electron_diffusion_length", GaAs_p) + GaAs_n.hole_diffusion_length = carrier_constants("hole_diffusion_length", GaAs_n) + + wls = np.linspace(300, 950, 100)*1e-9 + options = State() + options.wavelength = wls + options.T = 300 + options.light_iv = True + options.da_mode = 'green' + options.optics_method = 'BL' + + mesh = np.linspace(0, 2500, 1000) * 1e-9 + + pn_junction = Junction([Layer(si('500nm'), GaAs_p, role='emitter'), + Layer(si('2000nm'), GaAs_n, role='base')], + kind='DA', R_shunt=0.1, mesh=mesh, + sn=1e5, sp=1e5) + + solar_cell_solver(SolarCell([pn_junction]), 'optics', options) + + qe_sesame(pn_junction, options) + + sesame_EQE = pn_junction.eqe(wls) + sesame_IQE = pn_junction.iqe(wls) + + qe_depletion(pn_junction, options) + depletion_EQE = pn_junction.eqe(wls) + depletion_IQE = pn_junction.iqe(wls) + + assert sesame_EQE == approx(depletion_EQE, rel=0.1) + assert sesame_IQE == approx(depletion_IQE, rel=0.1) + + +def test_mesh_generation(): + from solcore.sesame_drift_diffusion.process_structure import make_mesh + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.state import State + + GaAs_n = material('GaAs')() + GaAs_p = material('GaAs')() + + options = State(minimum_spacing=1e-9, maximum_spacing=1e-9) + + pn_junction = Junction([Layer(si('500nm'), GaAs_n), + Layer(si('2000nm'), GaAs_p)]) + + layer_width = [layer.width*1e2 for layer in pn_junction] + + make_mesh(pn_junction, layer_width, options, [500e-7]) + + assert pn_junction.mesh == approx(np.arange(0, 2500.01, 1)*1e-9) + + +def test_get_srv_np(np_junction): + from solcore.sesame_drift_diffusion.process_structure import get_srv, process_structure + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.state import State + + options = State(T=300) + + GaAs_n = material('GaAs')(T=300, Nd=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_p = material('GaAs')(T=300, Na=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + + junction = Junction([Layer(si('200nm'), GaAs_n, role='emitter'), Layer(si('2000nm'), GaAs_p, role='base')]) + + junction.sn = 5 + junction.sp = 8 + + process_structure(junction, options) + + Sfront_e, Sfront_h, Sback_e, Sback_h = get_srv(junction) + + assert Sfront_e == junction.sn*100 + assert Sfront_h == junction.sn*100 + + assert Sback_e == junction.sp*100 + assert Sback_e == junction.sp*100 + + junction.sn_e = 2 + junction.sn_h = 3 + junction.sp_e = 4 + junction.sp_h = 5 + + Sfront_e, Sfront_h, Sback_e, Sback_h = get_srv(junction) + + assert Sfront_e == junction.sn_e*100 + assert Sfront_h == junction.sn_h*100 + + assert Sback_e == junction.sp_e*100 + assert Sback_h == junction.sp_h*100 + +def test_get_srv_pn(np_junction): + from solcore.sesame_drift_diffusion.process_structure import get_srv, process_structure + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.state import State + + options = State(T=300) + + GaAs_p = material('GaAs')(T=300, Na=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_n = material('GaAs')(T=300, Nd=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + + junction = Junction([Layer(si('200nm'), GaAs_p, role='emitter'), Layer(si('2000nm'), GaAs_n, role='base')]) + + junction.sn = 5 + junction.sp = 8 + + process_structure(junction, options) + + Sfront_e, Sfront_h, Sback_e, Sback_h = get_srv(junction) + + assert Sfront_e == junction.sp*100 + assert Sfront_h == junction.sp*100 + + assert Sback_e == junction.sn*100 + assert Sback_e == junction.sn*100 + + junction.sn_e = 2 + junction.sn_h = 3 + junction.sp_e = 4 + junction.sp_h = 5 + + Sfront_e, Sfront_h, Sback_e, Sback_h = get_srv(junction) + + assert Sfront_e == junction.sp_e*100 + assert Sfront_h == junction.sp_h*100 + + assert Sback_e == junction.sn_e*100 + assert Sback_h == junction.sn_h*100 + +def test_voltages_np(): + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.solar_cell_solver import solar_cell_solver, SolarCell + from solcore.sesame_drift_diffusion import iv_sesame + from solcore.state import State + from solcore.light_source import LightSource + + GaAs_n = material('GaAs')(T=300, Nd=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_p = material('GaAs')(T=300, Na=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + + options = State() + options.wavelength = np.linspace(300, 950, 100)*1e-9 + options.T = 300 + options.light_iv = True + options.light_source = LightSource(source_type='standard', version='AM1.5g', x=options.wavelength, output_units='photon_flux_per_m') + options.optics_method = 'TMM' + + np_junction = Junction([Layer(si('200nm'), GaAs_n, role='emitter'), Layer(si('2000nm'), GaAs_p, role='base')], + kind='sesame_PDD') + + solar_cell_solver(SolarCell([np_junction]), 'optics', options) + + voltage_end = [0, 1.3, 1.3] + voltage_points = [20, 41, 40] + + interp_voltages = np.linspace(-1, 0, 10) + + interp_results = np.zeros((len(voltage_end), len(interp_voltages))) + + for i, V_end in enumerate(voltage_end): + options.voltages = np.linspace(-1.3, V_end, voltage_points[i]) + options.internal_voltages = np.linspace(-1.3, V_end, voltage_points[i]) + + iv_sesame(np_junction, options) + + interp_results[i] = np_junction.iv(interp_voltages) + + assert interp_results[0] == approx(interp_results[1], rel=0.02) + assert interp_results[1] == approx(interp_results[2], rel=0.02) + + +def test_voltages_pn(): + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.solar_cell_solver import solar_cell_solver, SolarCell + from solcore.sesame_drift_diffusion import iv_sesame + from solcore.state import State + from solcore.light_source import LightSource + + GaAs_p = material('GaAs')(T=300, Na=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_n = material('GaAs')(T=300, Nd=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + + options = State() + options.wavelength = np.linspace(300, 950, 100)*1e-9 + options.T = 300 + options.light_iv = True + options.light_source = LightSource(source_type='standard', version='AM1.5g', x=options.wavelength, output_units='photon_flux_per_m') + options.optics_method = 'TMM' + + pn_junction = Junction([Layer(si('200nm'), GaAs_p, role='emitter'), Layer(si('2000nm'), GaAs_n, role='base')], + kind='sesame_PDD') + + solar_cell_solver(SolarCell([pn_junction]), 'optics', options) + + voltage_start = [0, -1.3, -1.3] + voltage_points = [20, 41, 40] + + interp_voltages = np.linspace(0, 1, 10) + + interp_results = np.zeros((len(voltage_start), len(interp_voltages))) + + for i, V_start in enumerate(voltage_start): + options.voltages = np.linspace(V_start, 1.3, voltage_points[i]) + options.internal_voltages = np.linspace(V_start, 1.3, voltage_points[i]) + + iv_sesame(pn_junction, options) + + interp_results[i] = pn_junction.iv(interp_voltages) + + assert interp_results[0] == approx(interp_results[1], rel=0.02) + assert interp_results[1] == approx(interp_results[2], rel=0.02) + + +def test_reverse_bias(): + from solcore import material, si + from solcore.structure import Junction, Layer + from solcore.solar_cell_solver import solar_cell_solver, SolarCell + from solcore.sesame_drift_diffusion import iv_sesame + from solcore.state import State + from solcore.light_source import LightSource + + GaAs_n = material('GaAs')(T=300, Nd=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + GaAs_p = material('GaAs')(T=300, Na=1e24, + hole_minority_lifetime=1e-6, electron_minority_lifetime=1e-6) + + options = State() + options.wavelength = np.linspace(300, 950, 100)*1e-9 + options.T = 300 + options.light_iv = True + options.light_source = LightSource(source_type='standard', version='AM1.5g', x=options.wavelength, output_units='photon_flux_per_m') + options.optics_method = 'TMM' + + np_junction = Junction([Layer(si('200nm'), GaAs_n, role='emitter'), Layer(si('2000nm'), GaAs_p, role='base')], + kind='sesame_PDD') + + solar_cell_solver(SolarCell([np_junction]), 'optics', options) + options.voltages = np.linspace(0, 1.3, 10) + options.internal_voltages = options.voltages + + iv_sesame(np_junction, options) + + # this is reverse bias for an n-p junction, so check that all currents are > 0: + + assert np.all(np_junction.current > 0) + assert np.all(np_junction.iv(options.voltages) > 0) + +