diff --git a/.github/workflows/build_deploy_wheels.yml b/.github/workflows/build_deploy_wheels.yml index 30fb4bd9..5655eee5 100644 --- a/.github/workflows/build_deploy_wheels.yml +++ b/.github/workflows/build_deploy_wheels.yml @@ -5,6 +5,7 @@ on: branches: - main - develop + tags: - "**" @@ -28,7 +29,8 @@ jobs: - [ macos-11, macosx, arm64] # cross compiled - [ windows-2019, win, AMD64 ] - python: [[ "cp37", "3.7" ], [ "cp38", "3.8" ], [ "cp39", "3.9" ], [ "cp310", "3.10" ], [ "cp311", "3.11" ]] + python: [[ "cp37", "3.7" ], [ "cp38", "3.8" ], [ "cp39", "3.9" ], + [ "cp310", "3.10" ], [ "cp311", "3.11" ], ["cp312", "3.12"]] name: Build wheel for ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }} ${{ matrix.buildplat[2] }} @@ -48,17 +50,18 @@ jobs: echo "c:\rtools40\ucrt64\bin;" >> $env:GITHUB_PATH - name: Install cibuildwheel - run: python -m pip install cibuildwheel==2.11.4 + run: python -m pip install cibuildwheel==2.16.2 - name: Build Solcore if: >- ( ! contains(matrix.buildplat[2], 'arm64' ) ) - uses: pypa/cibuildwheel@v2.11.4 + uses: pypa/cibuildwheel@v2.16.2 env: CIBW_BUILD: ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }}* CIBW_ARCHS: ${{ matrix.buildplat[2] }} CIBW_ENVIRONMENT_PASS_LINUX: RUNNER_OS CIBW_BEFORE_BUILD_MACOS: brew reinstall gfortran + CIBW_BEFORE_BUILD_LINUX: python -m pip install numpy --config-settings=setup-args="-Dallow-noblas=true" - name: Set extra env for arm64 if: >- @@ -69,7 +72,7 @@ jobs: - name: Cross-build Solcore for arm64 if: ${{ (matrix.python[0] != 'cp37') && ( contains(matrix.buildplat[2], 'arm64') )}} # image not present for python3.7 - uses: pypa/cibuildwheel@v2.11.4 + uses: pypa/cibuildwheel@v2.16.2 env: CIBW_BUILD: ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }}* CIBW_ARCHS: ${{ matrix.buildplat[2] }} diff --git a/.github/workflows/test_unit_and_examples.yml b/.github/workflows/test_unit_and_examples.yml index 079c513a..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.7", "3.8", "3.9", "3.10" ] + 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: @@ -98,7 +99,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [ "3.7", "3.8", "3.9", "3.10" ] + python-version: ["3.9", "3.10", "3.11", "3.12"] steps: - name: Checkout @@ -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 81645c89..1528fb8b 100755 --- a/docs/source/Solvers/DDsolver.rst +++ b/docs/source/Solvers/DDsolver.rst @@ -1,10 +1,22 @@ -Poisson - Drift-Diffusion solver (PDD) -====================================== +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`. 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 new file mode 100644 index 00000000..f39682b4 --- /dev/null +++ b/docs/source/Solvers/SesameDDsolver.rst @@ -0,0 +1,91 @@ +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 new file mode 100644 index 00000000..538e3020 --- /dev/null +++ b/docs/source/Solvers/sign_conventions.rst @@ -0,0 +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 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_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 42424b48..d1b0dab4 100755 --- a/examples/MJ_solar_cell_PDD_solver.py +++ b/examples/MJ_solar_cell_PDD_solver.py @@ -1,94 +1,153 @@ -import matplotlib.pyplot as plt import numpy as np +import matplotlib.pyplot as plt -from solcore.solar_cell import SolarCell, default_GaAs -from solcore.structure import Layer, Junction -from solcore import si -from solcore import material -from solcore.light_source import LightSource +from solcore import siUnits, material, si +from solcore.solar_cell import SolarCell +from solcore.structure import Junction, Layer from solcore.solar_cell_solver import solar_cell_solver +from solcore.light_source import LightSource -# TODO The purpose of this example needs to be clarified. It is messy - -T = 298 - -substrate = material('GaAs')(T=T) - - -def AlGaAs(T): - # We create the other materials we need for the device - window = material('AlGaAs')(T=T, Na=5e24, Al=0.8) - p_AlGaAs = material('AlGaAs')(T=T, Na=1e24, Al=0.4) - n_AlGaAs = material('AlGaAs')(T=T, Nd=8e22, Al=0.4) - bsf = material('AlGaAs')(T=T, Nd=2e24, Al=0.6) - - output = 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') - - return output - - -my_solar_cell = SolarCell([default_GaAs(T)], T=T, R_series=0, substrate=substrate) - -Vin = np.linspace(-2, 2.61, 201) -V = np.linspace(0, 2.6, 300) -wl = np.linspace(350, 1000, 301) * 1e-9 -light_source = LightSource(source_type='standard', version='AM1.5g', x=wl, output_units='photon_flux_per_m', - concentration=1) - -if __name__ == '__main__': - # We calculate the IV curve under illumination - solar_cell_solver(my_solar_cell, 'equilibrium') - # , - # user_options={'T_ambient': T, 'db_mode': 'boltzmann', 'voltages': V, 'light_iv': True, - # 'wavelength': wl, 'optics_method': 'BL', 'mpp': True, 'internal_voltages': Vin, - # 'light_source': light_source}) - - # We can plot the electron and hole densities in equilibrium and at short circuit, both calculated automatically - # before calculating the IV curve - plt.figure(1) - for j in my_solar_cell.junction_indices: - zz = my_solar_cell[j].equilibrium_data.Properties['x'] + my_solar_cell[j].offset - n = my_solar_cell[j].equilibrium_data.Properties['Nd'] - p = my_solar_cell[j].equilibrium_data.Properties['Na'] - plt.semilogy(zz, n, 'b') - plt.semilogy(zz, p, 'r') - - zz = my_solar_cell[j].equilibrium_data.Bandstructure['x'] + my_solar_cell[j].offset - n = my_solar_cell[j].equilibrium_data.Bandstructure['n'] - p = my_solar_cell[j].equilibrium_data.Bandstructure['p'] - plt.semilogy(zz, n, 'b--') - plt.semilogy(zz, p, 'r--') - - plt.xlabel('Position (m)') - plt.ylabel('Carrier density (m$^{-3}$)') - # - # # And the IV curves of the individual junctions and of the MJ device - # plt.figure(2) - # plt.plot(V, abs(my_solar_cell[0].iv(V)), 'b', label='AlGaAs') - # plt.plot(V, abs(my_solar_cell[1].iv(V)), 'r', label='GaAs') - # plt.plot(my_solar_cell.iv.IV[0], abs(my_solar_cell.iv.IV[1]), 'k', label='MJ') - # - # plt.legend() - # plt.xlim(0, 2.6) - # plt.ylim(0, 200) - # plt.xlabel('Bias (V)') - # plt.ylabel('Current (A/m$^2}$)') - # - # # Now we calculate the quantum efficiency - # solar_cell_solver(my_solar_cell, 'qe', - # user_options={'T_ambient': T, 'db_mode': 'boltzmann', 'voltages': V, 'light_iv': True, - # 'wavelength': wl, 'optics_method': 'BL', 'mpp': True, 'internal_voltages': Vin, - # 'light_source': light_source}) - # - # plt.figure(3) - # - # plt.plot(wl * 1e9, my_solar_cell[0].eqe(wl), 'b', label='AlGaAs') - # plt.plot(wl * 1e9, my_solar_cell[1].eqe(wl), 'r', label='GaAs') - # - # plt.xlabel('Wavelength (nm)') - # plt.ylabel('EQE') - - plt.show() +# Define materials for the anti-reflection coating: +MgF2 = material("MgF2")() +ZnS = material("ZnScub")() + +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 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(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) + +# MID CELL - GaAs +GaAs = material("GaAs") + +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"), + 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 = '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 + + [ + Junction([ + Layer(si("25nm"), material=window_material, role='window'), + Layer(si("100nm"), material=top_cell_n_material, role='emitter'), + 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'), + ], sn=1e4, sp=1e4, kind='PDD'), + Junction([Layer(si("400nm"), material=bot_cell_n_material, role='emitter'), + Layer(si("100um"), material=bot_cell_p_material, role='base'), + ], sn=1e4, sp=1e4, kind='PDD') + ], + shading=0.02, cell_area=1 * 1 / 1e4) + +# Choose wavelength range (in m): +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, + '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(-2.8, 0, 50) +internal_voltages = np.linspace(-3, 1.7, 120) + +# 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(-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(-10, 35) +plt.xlim(0, 3) +plt.ylabel('Current (mA/cm$^2$)') +plt.xlabel('Voltage (V)') +plt.show() \ No newline at end of file diff --git a/examples/MJ_solar_cell_using_DA.py b/examples/MJ_solar_cell_using_DA.py index 948e0cdd..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,51 +80,66 @@ # 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') -V = np.linspace(0, 3, 300) + +# 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) # 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, - # '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(V, solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell') -plt.plot(V, -solar_cell[0 + len(ARC_layers)].iv(V)/10, 'b', label='InGaP sub-cell') -plt.plot(V, -solar_cell[1 + len(ARC_layers)].iv(V)/10, 'g', label='GaAs sub-cell') -plt.plot(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.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} %') diff --git a/examples/Optical_constants_RAT_of_ARC.py b/examples/Optical_constants_RAT_of_ARC.py index 9a605fa3..51c13c24 100755 --- a/examples/Optical_constants_RAT_of_ARC.py +++ b/examples/Optical_constants_RAT_of_ARC.py @@ -1,68 +1,62 @@ -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 +wl = np.linspace(300, 1900, 1000) -E_eV = np.linspace(0.65, 4, 1000) +MgF2 = material("MgF2")() +HfO2 = material("HfO2")() +ZnS = material("ZnScub")() +AlInP = material("AlInP")(Al=0.52) +GaInP = material("GaInP")(In=0.49) - -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!") - - # 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) - - return (energy, n_interp, k_interp) - - -# Load nk data using nk_convert function... -alinp_nk = nk_convert("data/AlInP_nk.csv", energy=E_eV) -gainp_nk = nk_convert("data/GaInP_nk.csv", energy=E_eV) -mgf_nk = nk_convert("data/MgF_nk.csv", energy=E_eV) -sic_nk = nk_convert("data/SiCx_nk.csv", energy=E_eV) -zns_nk = nk_convert("data/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/examples/Si_cell.py b/examples/Si_cell.py new file mode 100644 index 00000000..b43663f6 --- /dev/null +++ b/examples/Si_cell.py @@ -0,0 +1,197 @@ +from solcore.solar_cell import SolarCell, Junction, Layer +from solcore.state import State +from solcore.solar_cell_solver import solar_cell_solver +import numpy as np +from solcore.light_source import LightSource +from scipy.special import erfc + +import matplotlib.pyplot as plt + +from solcore import material, si + +d_bulk = 100e-6 + +Air = material("Air")() +MgF2 = material("MgF2")() +Ag = material("Ag")() + +Si_pn = material("Si")(electron_mobility=si("1e4cm2"), hole_mobility=si("1e3cm2"), + electron_minority_lifetime=0.001, hole_minority_lifetime=0.001) + +wavelengths = np.linspace(280, 1200, 100)*1e-9 + +TCO = material('ITO2')() + +front_materials = [Layer(80e-9, MgF2), Layer(55e-9, TCO)] + +back_materials = [Layer(55e-9, TCO), + Layer(120e-9, MgF2)] + +options = State() +options.wavelength = wavelengths +options.optics_method = 'TMM' +options.light_iv = True +options.T = 298 + +options.light_source = LightSource(source_type="standard", + version="AM1.5g", x=options.wavelength) + +options.voltages = np.linspace(0, 0.8, 40) +options.internal_voltages = options.voltages +options.mpp = True +options.no_back_reflection = False +options.position = 10e-9 + +nD = si("1e20cm-3") +nA = si("1e19cm-3") +bulk_doping = si("5e15cm-3") # n type bulk + +# rear junction (n-type) +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 + +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', + **surface_recomb, + ) + ] + +Si_cell = SolarCell(front_materials + + Si_junction + + back_materials, + shading=0.02, + substrate=Ag, + ) + +solar_cell_solver(Si_cell, 'iv', options) + +solar_cell_solver(Si_cell, 'qe', options) + +result_stack = np.vstack([Si_cell.reflected, [layer.layer_absorption for layer in Si_cell]]) + +fig, (ax, ax2) = plt.subplots(1, 2, figsize=(10, 3.5)) +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), '-k', linewidth=2, + label='EQE') + +ax.set_xlim(280, 1200) +ax.set_ylim(0, 100) +ax.set_xlabel("Wavelength (nm)") +ax.set_ylabel("R / A / EQE (%)") +ax.set_title('a) EQE and cell optics', loc='left') +ax.legend() +# plt.show() + +jsc = Si_cell.iv.Isc / 10 + +ax2.plot(Si_cell.iv['IV'][0], Si_cell.iv['IV'][1] / 10, '-', label='IV', + linewidth=2, color='k') + +ax2.set_ylim(0, 1.03 * jsc) +ax2.set_xlim(np.min(options.voltages), np.max(options.voltages)) +ax2.set_xlabel('Voltage (V)') +ax2.set_ylabel('Current density (mA/cm$^2$)') +ax2.set_title('b) IV characteristics and power output', loc='left') + +ax3 = ax2.twinx() +ax3.plot(options.voltages, Si_cell.iv['IV'][0] * Si_cell.iv['IV'][1], + '-r', label='Power', linewidth=2) +ax3.set_ylabel('Power density (W m$^{-2}$)') +ax3.set_ylim(0, 1.03 * jsc * 10) + +ax3.spines['right'].set_color('r') +ax3.yaxis.label.set_color('r') +ax3.tick_params(axis='y', colors='r') + +ax2.set_axisbelow(True) +ax3.set_axisbelow(True) + +ax2.text(0.02, 0.9 * jsc, r'$J_{SC}$', zorder=5) +ax2.text(0.02, 0.8 * jsc, r'$V_{OC}$') +ax2.text(0.02, 0.7 * jsc, 'FF') +ax2.text(0.02, 0.6 * jsc, r'$\eta$') +ax2.text(0.02, 0.5 * jsc, r'$J_{MPP}$') +ax2.text(0.02, 0.4 * jsc, r'$V_{MPP}$') + +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.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) + +ax3.grid(False) +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/pyproject.toml b/pyproject.toml index 705a2e70..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,6 +39,7 @@ dependencies = [ "pyyaml", "yabox", "joblib", + "solsesame", ] dynamic = ["version"] @@ -63,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/IV.py b/solcore/analytic_solar_cells/IV.py index c55ab0ae..24843d5f 100755 --- a/solcore/analytic_solar_cells/IV.py +++ b/solcore/analytic_solar_cells/IV.py @@ -11,7 +11,7 @@ def iv_multijunction(solar_cell, options): points. If photocurrent is not provided, the resultant IV characteristics are purely recombination currents, otherwise light IVs are returned. - In the end, the soalr_cell object is updated with an "iv" attribute containing a dictionary with: + In the end, the solar_cell object is updated with an "iv" attribute containing a dictionary with: "IV": (V, I) Calculated IV characteristics "junction IV": [(V junc 1, I junc 1), (V junc 2, I junc 2), ...] "Rseries IV": (V, I) Calculated IV characteristics of the series resistance @@ -103,8 +103,8 @@ def iv_multijunction(solar_cell, options): try: Isc = np.interp([0], VV, II)[0] Voc = np.interp([0], -II, VV)[0] - Pmpp = np.max(PP) - idx = np.argmax(PP) + Pmpp = np.nanmax(PP) + idx = np.nanargmax(PP) Vmpp = VV[idx] Impp = II[idx] FF = Pmpp / (Isc * Voc) diff --git a/solcore/analytic_solar_cells/depletion_approximation.py b/solcore/analytic_solar_cells/depletion_approximation.py index 22ee67c6..37bbef4f 100755 --- a/solcore/analytic_solar_cells/depletion_approximation.py +++ b/solcore/analytic_solar_cells/depletion_approximation.py @@ -160,11 +160,29 @@ def iv_depletion(junction, options): "J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003).", ) - junction.voltage = options.internal_voltages T = options.T kbT = kb * T id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(junction) + + if pn_or_np == "pn": + junction.voltage = options.internal_voltages + J_sign = 1 + + else: + junction.voltage = -options.internal_voltages + J_sign = -1 + + if np.all(options.internal_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) + + if "voltages" in options.keys(): + 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) + + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es_n, es_p = identify_parameters( junction, T, pRegion, nRegion, iRegion ) @@ -307,7 +325,10 @@ def iv_depletion(junction, options): J_sc_scr = q * get_J_sc_SCR(xa, xb, g, wl, ph) # And, finally, we output the currents - junction.current = ( + junction.voltage = J_sign*junction.voltage # this flips the voltage sign back for an n-p junction, or leaves it + # unchanged for a p-n junction + + junction.current = J_sign*( Jrec + JnDark + JpDark + V / R_shunt - J_sc_top - J_sc_bot - J_sc_scr ) junction.iv = interp1d( @@ -320,12 +341,12 @@ def iv_depletion(junction, options): ) junction.region_currents = State( { - "Jn_dif": JnDark, - "Jp_dif": JpDark, - "Jscr_srh": Jrec, - "J_sc_top": J_sc_top, - "J_sc_bot": J_sc_bot, - "J_sc_scr": J_sc_scr, + "Jn_dif": J_sign*JnDark, + "Jp_dif": J_sign*JpDark, + "Jscr_srh": J_sign*Jrec, + "J_sc_top": J_sign*J_sc_top, + "J_sc_bot": J_sign*J_sc_bot, + "J_sc_scr": J_sign*J_sc_scr, } ) @@ -343,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/Levinshtein/GroupIV.txt b/solcore/material_data/Levinshtein/GroupIV.txt index 82e26c1f..79c02d14 100755 --- a/solcore/material_data/Levinshtein/GroupIV.txt +++ b/solcore/material_data/Levinshtein/GroupIV.txt @@ -25,6 +25,8 @@ electron_affinity=4.0 eV electron_auger_recombination=1e-42 hole_auger_recombination=1e-42 radiative_recombination=6.4e-20 +Nc=1e25 +Nv=5e24 [Si] lattice_constant=5.431 Angstrom @@ -52,4 +54,6 @@ relative_permittivity=11.7 electron_affinity=4.05 eV electron_auger_recombination=1.1e-42 hole_auger_recombination=3e-43 -radiative_recombination=1.1e-20 \ No newline at end of file +radiative_recombination=1.1e-20 +Nc=3e25 +Nv=1e25 \ No newline at end of file diff --git a/solcore/material_data/Vurgaftman/ternaries.txt b/solcore/material_data/Vurgaftman/ternaries.txt index 1760a429..77dd650d 100755 --- a/solcore/material_data/Vurgaftman/ternaries.txt +++ b/solcore/material_data/Vurgaftman/ternaries.txt @@ -90,7 +90,7 @@ d=0 eV parent0=GaAs parent1=GaN x=N -Eg_Gamma=120.4-100.*x eV +Eg_Gamma=120.4-100.*N eV spin_orbit_splitting=0 eV [AlGaP] 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 48e47c29..1f3799ac 100755 --- a/solcore/material_system/material_system.py +++ b/solcore/material_system/material_system.py @@ -326,22 +326,31 @@ def __getattr__(self, attrname): # only used for unknown attributes. except: return calculate_mobility(self.material_string, True, self.Na, self.main_fraction) if attrname == "Nc": + try: + return ParameterSystem().get_parameter(self.material_string, attrname) # return 2 * (2 * pi * self.eff_mass_electron_Gamma * m0 * kb * self.T / h ** 2) ** 1.5 # Changed by Diego 22/10/18: some materials dont have _Gamma but all have the normal electron effective mass - return 2 * (2 * pi * self.eff_mass_electron * m0 * kb * self.T / h ** 2) ** 1.5 + except: + return 2 * (2 * pi * self.eff_mass_electron * m0 * kb * self.T / h ** 2) ** 1.5 if attrname == "Nv": + try: + return ParameterSystem().get_parameter(self.material_string, attrname) + except: # Strictly speaking, this is valid only for III-V, zinc-blend semiconductors - mhh = self.eff_mass_hh_z - mlh = self.eff_mass_lh_z - Nvhh = 2 * (2 * pi * mhh * m0 * kb * self.T / h ** 2) ** 1.5 - Nvlh = 2 * (2 * pi * mlh * m0 * kb * self.T / h ** 2) ** 1.5 - return Nvhh + Nvlh + mhh = self.eff_mass_hh_z + mlh = self.eff_mass_lh_z + Nvhh = 2 * (2 * pi * mhh * m0 * kb * self.T / h ** 2) ** 1.5 + Nvlh = 2 * (2 * pi * mlh * m0 * kb * self.T / h ** 2) ** 1.5 + return Nvhh + Nvlh if attrname == "ni": return np.sqrt(self.Nc * self.Nv * np.exp(-self.band_gap / (kb * self.T))) if attrname == "radiative_recombination": - inter = lambda E: self.n(E) ** 2 * self.alphaE(E) * np.exp(-E / (kb * self.T)) * E ** 2 - upper = self.band_gap + 10 * kb * self.T - return 1.0 / self.ni ** 2 * 2 * pi / (h ** 3 * c ** 2) * quad(inter, 0, upper)[0] + try: + return ParameterSystem().get_parameter(self.material_string, attrname) + except: + inter = lambda E: self.n(E) ** 2 * self.alphaE(E) * np.exp(-E / (kb * self.T)) * E ** 2 + upper = self.band_gap + 10 * kb * self.T + 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 diff --git a/solcore/meson.build b/solcore/meson.build index 42e1a582..ac5bc589 100644 --- a/solcore/meson.build +++ b/solcore/meson.build @@ -9,6 +9,7 @@ install_subdir('optics', install_dir: py.get_install_dir() / 'solcore') install_subdir('optimization', install_dir: py.get_install_dir() / 'solcore') install_subdir('parameter_system', install_dir: py.get_install_dir() / 'solcore') subdir('poisson_drift_diffusion') +install_subdir('sesame_drift_diffusion', install_dir: py.get_install_dir() / 'solcore') install_subdir('quantum_mechanics', install_dir: py.get_install_dir() / 'solcore') install_subdir('spice', install_dir: py.get_install_dir() / 'solcore') install_subdir('units_system', install_dir: py.get_install_dir() / 'solcore') 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/DDmodel-current.f95 b/solcore/poisson_drift_diffusion/DDmodel-current.f95 old mode 100755 new mode 100644 diff --git a/solcore/poisson_drift_diffusion/DeviceStructure.py b/solcore/poisson_drift_diffusion/DeviceStructure.py old mode 100755 new mode 100644 diff --git a/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py b/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py old mode 100755 new mode 100644 index 12fe9de1..3667559a --- a/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py +++ b/solcore/poisson_drift_diffusion/DriftDiffusionUtilities.py @@ -4,15 +4,15 @@ from numpy.typing import NDArray from scipy.interpolate import interp1d -from .. import asUnit, constants -from ..light_source import LightSource -from ..registries import ( +from solcore import asUnit, constants +from solcore.light_source import LightSource +from solcore.registries import ( register_equilibrium_solver, register_iv_solver, register_short_circuit_solver, ) -from ..state import State -from ..structure import Junction +from solcore.state import State +from solcore.structure import Junction from .DeviceStructure import CalculateAbsorptionProfile, CreateDeviceStructure try: @@ -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/poisson_drift_diffusion/QWunit.py b/solcore/poisson_drift_diffusion/QWunit.py old mode 100755 new mode 100644 diff --git a/solcore/poisson_drift_diffusion/__init__.py b/solcore/poisson_drift_diffusion/__init__.py old mode 100755 new mode 100644 index 18134fb8..5e7fc792 --- a/solcore/poisson_drift_diffusion/__init__.py +++ b/solcore/poisson_drift_diffusion/__init__.py @@ -9,5 +9,4 @@ except Exception as err: print('WARNING: The Poisson - Drift-Diffusion solver will not be available because ' 'the ddModel fortran library could not be imported.') - print(err) - + print(err) \ No newline at end of file diff --git a/solcore/poisson_drift_diffusion/meson.build b/solcore/poisson_drift_diffusion/meson.build index 4973accf..da6df704 100644 --- a/solcore/poisson_drift_diffusion/meson.build +++ b/solcore/poisson_drift_diffusion/meson.build @@ -68,7 +68,7 @@ if get_option('with_pdd') input : ['DDmodel-current.f95'], output : ['ddModelmodule.c', 'ddModel-f2pywrappers2.f90'], command : [py, '-m', 'numpy.f2py', '-m', 'ddModel', '@INPUT@', '--lower', - '--build-dir', 'solcore/poisson_drift_diffusion'] + '--build-dir', 'solcore/poisson_drift_diffusion/'] ) py.extension_module('ddModel', @@ -92,5 +92,5 @@ python_sources = [ py.install_sources( python_sources, - subdir: 'solcore' / 'poisson_drift_diffusion' + subdir: 'solcore' / 'poisson_drift_diffusion', ) diff --git a/solcore/sesame_drift_diffusion/__init__.py b/solcore/sesame_drift_diffusion/__init__.py new file mode 100644 index 00000000..c4b27a7f --- /dev/null +++ b/solcore/sesame_drift_diffusion/__init__.py @@ -0,0 +1 @@ +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 new file mode 100644 index 00000000..b6add5f7 --- /dev/null +++ b/solcore/sesame_drift_diffusion/solve_pdd.py @@ -0,0 +1,420 @@ +from __future__ import annotations + +from sesame import Builder, IVcurve, Analyzer +import numpy as np +from scipy.optimize import root +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 solcore.registries import register_iv_solver, register_equilibrium_solver + +@register_equilibrium_solver("sesame_PDD") +def equilibrium(junction: Junction, + **kwargs): + """ + Solve at equilibrium (no illumination, no applied voltage) using the Sesame solver. + + :param junction: a Junction object + :param options: a State object containing options for the solver + """ + + # TODO: pass output from 'result' parameter to the user, similar to the Fortran PDD solver, so that + # band structures etc. can be plotted. + + options = State(kwargs) # needs to be passed as kwargs to be compatible with Fortran equilibrium solver + # in registries + + if not hasattr(junction, "sesame_sys"): + process_structure(junction, options) + + j, result = IVcurve(junction.sesame_sys, [0]) + + if np.any(np.isnan(j)): + warnings.warn( + "Current calculation did not converge at all voltages", UserWarning + ) + + j = j * junction.sesame_sys.scaling.current * 1e4 # cm-2 -> m-2 + + junction.sesame_output = result + + +@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. + + :param junction: a Junction object + :param options: a State object containing options for the solver + """ + + # 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) + + 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 + ) + + 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? + junction.sesame_sys.generation(g_vs_z) + + voltages = options.internal_voltages + + R_shunt = min(junction.R_shunt, 1e14) if hasattr(junction, "R_shunt") else 1e14 + + 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: + # 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, + ) + + 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): + # positive and negative voltages + + negative_voltages = voltages_for_solve[voltages_for_solve <= 0] + positive_voltages = voltages_for_solve[voltages_for_solve >= 0] + + negative_voltages_order = np.argsort(negative_voltages)[::-1] + positive_voltages_order = np.argsort(positive_voltages) + + negative_voltages = negative_voltages[negative_voltages_order] + positive_voltages = positive_voltages[positive_voltages_order] + + 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_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) + ) + + else: + j = np.concatenate((j_negative, j_positive)) + 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 + else: + # negative voltages only + voltage_order = np.argsort(voltages_for_solve)[::-1] + final_voltages = voltages_for_solve[voltage_order] + 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 = 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: + result_voltage = -final_voltages + sort_result = np.argsort(result_voltage) + j = j[sort_result] + result = {key: result[key][sort_result, :] for key in result.keys()} + result_voltage = result_voltage[sort_result] + + else: + result_voltage = final_voltages + + 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] + + 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]] + + 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[non_nans], + shunted_current[non_nans], + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(first_non_nan, last_non_nan), + ) + + 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. + + :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) + + # if not options.parallel: + # n_jobs = 1 + # + # else: + # n_jobs = options.n_jobs if hasattr(options, "n_jobs") else -1 + # n_jobs = 1 + # parallel implementation does not work due to pickling error. Unsure of cause. + + wls = options.wavelength + + 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) + + 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(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 + + 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 + + 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)) + + if i1 == wl_solve[0]: + guess = None + + else: + guess = {key: result[key][0, :] for key in result.keys()} + + j, result = IVcurve(junction.sesame_sys, voltages, guess=guess) + + eqe[i1] = np.abs(j) / (q * flux) + + if np.any(np.isnan(eqe)): + warnings.warn( + "EQE calculation did not converge at all wavelengths", UserWarning + ) + + warnings.resetwarnings() + + eqe = eqe * junction.sesame_sys.scaling.current + 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 + + junction.iqe = interp1d(wls, iqe) + + junction.eqe = interp1d( + wls, + eqe, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(eqe[0], eqe[-1]), + ) + + 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/solcore/solar_cell.py b/solcore/solar_cell.py index 4d6b29ca..4a9895b4 100755 --- a/solcore/solar_cell.py +++ b/solcore/solar_cell.py @@ -1,8 +1,6 @@ from solcore.structure import Layer, Junction, Structure, TunnelJunction from solcore import material, si -import numpy as np - def default_GaAs(T): # We create the other materials we need for the device @@ -124,34 +122,3 @@ def update_junction(self, junction, **kwargs): def __call__(self, i): return self[self.junction_indices[i]] - - -if __name__ == '__main__': - window = material('AlGaAs')(T=298, Na=1e24, Al=0.8) - stack = [Layer(width=si("50nm"), material=window), - default_GaAs(298)] - - # stack = [default_GaAs(298)] - - my_cell = SolarCell(layers=stack) - # - # my_cell.append_multiple([default_GaAs(298), default_GaAs(298), default_GaAs(298)]) - # print(my_cell) - - from solcore.poisson_drift_diffusion.DriftDiffusionUtilities import solve_pdd, default_photon_flux, \ - default_wavelengths - import matplotlib.pyplot as plt - - solve_pdd(my_cell, 'QE', vfin=1.2, vstep=0.05, light=True) - QE = my_cell(0).qe - - # Finally, we plot the internal and external quantum efficiencies using the information stored in the output dictionaries - plt.plot(QE['QE']['wavelengths'] / 1e-9, QE['QE']['IQE'] * 100, label='IQE') - plt.plot(QE['QE']['wavelengths'] / 1e-9, QE['QE']['EQE'] * 100, label='EQE') - plt.plot(QE['QE']['wavelengths'] / 1e-9, my_cell.T * 100, label='T') - plt.ylim(0, 100) - plt.legend(loc='lower left') - plt.ylabel('QE (%)') - plt.xlabel('Wavelength (nm)') - - plt.show() diff --git a/solcore/solar_cell_solver.py b/solcore/solar_cell_solver.py index c6c72d8d..d929df50 100755 --- a/solcore/solar_cell_solver.py +++ b/solcore/solar_cell_solver.py @@ -4,7 +4,7 @@ import numpy as np -from . import analytic_solar_cells as ASC +from . import analytic_solar_cells as ASC, sesame_drift_diffusion as sesame_PDD, poisson_drift_diffusion as PDD from .light_source import LightSource from .optics import ( # noqa rcwa_options, @@ -25,16 +25,15 @@ from .structure import Junction, Layer, TunnelJunction try: - from . import poisson_drift_diffusion as PDD a = PDD.pdd_options except AttributeError: PDD.pdd_options = {} default_options = State() -pdd_options = PDD.pdd_options -asc_options = ASC.db_options - +# pdd_options = PDD.pdd_options +# asc_options = ASC.db_options +# sesame_options = sesame_PDD.sesame_options def merge_dicts(*dict_args): """ @@ -73,7 +72,7 @@ def merge_dicts(*dict_args): default_options.recalculate_absorption = False default_options = merge_dicts( - default_options, ASC.db_options, ASC.da_options, PDD.pdd_options, rcwa_options + default_options, ASC.db_options, ASC.da_options, PDD.pdd_options, rcwa_options, ) @@ -169,6 +168,8 @@ def solve_iv(solar_cell, options): if solar_cell[j].kind == "PDD": PDD.iv_pdd(solar_cell[j], **options) + elif solar_cell[j].kind == "sesame_PDD": + sesame_PDD.iv_sesame(solar_cell[j], options) elif solar_cell[j].kind == "DA": ASC.iv_depletion(solar_cell[j], options) elif solar_cell[j].kind == "2D": @@ -224,6 +225,8 @@ def solve_qe(solar_cell, options): for j in solar_cell.junction_indices: if solar_cell[j].kind == "PDD": PDD.qe_pdd(solar_cell[j], options) + elif solar_cell[j].kind == "sesame_PDD": + sesame_PDD.qe_sesame(solar_cell[j], options) elif solar_cell[j].kind == "DA": ASC.qe_depletion(solar_cell[j], options) elif solar_cell[j].kind == "2D": @@ -347,6 +350,18 @@ def prepare_solar_cell(solar_cell, options): process_position(solar_cell, options, layer_widths) + # check that the voltages are in the expected format (monotonically increasing): + if 'voltages' in options.keys(): + if np.any(np.diff(options.voltages) < 0): + raise ValueError( + "ERROR in 'solar_cell_solver' - The voltages must be monotonically increasing." + ) + + if 'internal_voltages' in options.keys(): + if np.any(np.diff(options.internal_voltages) < 0): + raise ValueError( + "ERROR in 'solar_cell_solver' - The internal_voltages must be monotonically increasing." + ) def process_position(solar_cell, options, layer_widths): """ diff --git a/tests/test__pdd.py b/tests/test__pdd.py index 6c4ba1a7..d3b04095 100755 --- a/tests/test__pdd.py +++ b/tests/test__pdd.py @@ -119,7 +119,7 @@ def test_calculate_iv(mocker): mock_dump_iv = mocker.patch(f"{root}.DumpIV") from solcore.structure import Junction - from solcore.poisson_drift_diffusion.DriftDiffusionUtilities import calculate_iv + from solcore.poisson_drift_diffusion import calculate_iv junction = Junction() vlimit = 1 @@ -151,7 +151,7 @@ def test_find_minimum_bandgap(): from solcore.structure import Layer, Junction from solcore import material from solcore.constants import q - from solcore.poisson_drift_diffusion.DriftDiffusionUtilities import ( + from solcore.poisson_drift_diffusion import ( find_minimum_bandgap, ) @@ -177,7 +177,7 @@ def kbT(T: float = 300) -> float: ], ) def test_find_voltage_limits(bandgap, p_on_n, vmax, vmin, exp_vmax, exp_vmin): - from solcore.poisson_drift_diffusion.DriftDiffusionUtilities import ( + from solcore.poisson_drift_diffusion import ( find_voltage_limits, ) @@ -187,7 +187,7 @@ def test_find_voltage_limits(bandgap, p_on_n, vmax, vmin, exp_vmax, exp_vmin): def test_consolidate_iv(): - from solcore.poisson_drift_diffusion.DriftDiffusionUtilities import ( + from solcore.poisson_drift_diffusion import ( consolidate_iv, ) diff --git a/tests/test_depletion_approximation.py b/tests/test_depletion_approximation.py index 5d8f566f..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)) ) ) @@ -1198,7 +1202,7 @@ def test_dark_iv_depletion_np(np_junction): options.light_iv = False T = options.T - test_junc[0].voltage = options.internal_voltages + test_junc[0].voltage = -options.internal_voltages id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers( test_junc[0] @@ -1242,9 +1246,9 @@ def test_dark_iv_depletion_np(np_junction): J_sc_bot = 0 J_sc_scr = 0 - current = Jrec + JnDark + JpDark + V / 1e14 - J_sc_top - J_sc_bot - J_sc_scr + current = -(Jrec + JnDark + JpDark + V / 1e14 - J_sc_top - J_sc_bot - J_sc_scr) iv = interp1d( - test_junc[0].voltage, + options.internal_voltages, current, kind="linear", bounds_error=False, 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 new file mode 100644 index 00000000..bc6f8a49 --- /dev/null +++ b/tests/test_sesame_pdd.py @@ -0,0 +1,514 @@ +import numpy as np +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.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) + Si_i = material('Si')(Na=0, Nd=0, 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) + + junction = Junction([Layer(si('200nm'), Si_n), + Layer(si('2000nm'), Si_i), + Layer(si('2000nm'), Si_p)]) + + options = State() + options.T = 300 + + process_structure(junction, options) + + sesame_obj = junction.sesame_sys + + assert len(np.unique(sesame_obj.rho)) == 3 + assert sesame_obj.rho[0] == 1e24*1e-6/sesame_obj.scaling.density + assert sesame_obj.rho[-1] == -2e24*1e-6/sesame_obj.scaling.density + + +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) + + 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 + 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) + } + + assert result == expected_dict + + 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_compare_DA_np(): + + 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.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, 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) + + 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) + options.internal_voltages = np.linspace(-1.3, 0.5, 30) + 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, 2200, 500)*1e-9 + + 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([np_junction]), 'optics', options) + + iv_depletion(np_junction, options) + depletion_current = np_junction.current + depletion_current_interp = np_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]) + + assert depletion_current_interp[-1] == approx(sesame_current_interp[-1], rel=0.05) + assert np.sign(depletion_current_interp[0]) == np.sign(sesame_current_interp[0]) + + assert np.all(sesame_current == sesame_current_interp) + + +def test_compare_DA_pn(): + + 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.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_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) + + options = State() + options.wavelength = np.linspace(300, 950, 100)*1e-9 + options.voltages = np.linspace(-0.5, 1.3, 30) + options.internal_voltages = np.linspace(-0.5, 1.3, 30) + 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, 2200, 500)*1e-9 + + 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) + + iv_depletion(pn_junction, options) + depletion_current = pn_junction.current + depletion_current_interp = pn_junction.iv(options.voltages) + + iv_sesame(pn_junction, options) + sesame_current = pn_junction.current + sesame_current_interp = pn_junction.iv(options.voltages) + + assert depletion_current[0] == approx(sesame_current[0], rel=0.05) + assert np.sign(depletion_current[-1]) == np.sign(sesame_current[-1]) + + assert depletion_current_interp[0] == approx(sesame_current_interp[0], rel=0.05) + 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) + + diff --git a/tests/test_structure.py b/tests/test_structure.py index 3c63eece..2bc5427a 100644 --- a/tests/test_structure.py +++ b/tests/test_structure.py @@ -1,7 +1,7 @@ from pytest import approx, fixture import random from solcore import material -from solcore.poisson_drift_diffusion.DeviceStructure import CreateDeviceStructure +from solcore.poisson_drift_diffusion import CreateDeviceStructure from solcore.structure import Junction, Layer, SolcoreMaterialToStr, Structure, ToSolcoreMaterial, TunnelJunction from solcore.structure import InLineComposition, ToLayer, ToStructure