Solve issue with boundary condition in depletion approximation

.devpy/cmds.py

@@ -1,6 +1,6 @@
 import click
 from devpy import util
-from devpy.cmds.util import get_site_packages, set_pythonpath, run
+from devpy.cmds.meson import _get_site_packages
 
 
 @click.command()
@@ -25,19 +25,19 @@ def install_dependencies(test_dep=False, doc_dep=False):
     config = util.get_config()
     default_dependencies = config["project"]['dependencies']
     print("Installing dependencies", default_dependencies)
-    run(
+    util.run(
         ["pip", "install"] + list(default_dependencies),
     )
     if test_dep:
         test_dependencies = config["project.optional-dependencies"]['test']
         print("Installing test-dependencies", test_dependencies)
-        run(
+        util.run(
             ["pip", "install"] + list(test_dependencies),
         )
     if doc_dep:
         doc_dependencies = config["project.optional-dependencies"]['doc']
         print("Installing doc-dependencies", doc_dependencies)
-        run(
+        util.run(
             ["pip", "install"] + list(doc_dependencies),
         )
 
@@ -60,9 +60,9 @@ def codecov(build_dir, codecov_args):
         codecov parameters
     """
 
-    site_path = get_site_packages(build_dir)
+    site_path = _get_site_packages()
 
-    run(
+    util.run(
         ["codecov"] + list(codecov_args),
         cwd=site_path,
-    )
+    )

.github/workflows/test_unit_and_examples.yml

@@ -53,7 +53,7 @@ jobs:
 
       - name: Install Python dependecies
         run: |
-          pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy
+          pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1
           python3 -m devpy install-dependencies -test-dep
 
       - name: Install S4
@@ -125,7 +125,7 @@ jobs:
 
       - name: Install Python dependecies
         run: |
-          pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy
+          pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1
           python3 -m devpy install-dependencies -test-dep
 
       - name: Install S4

docs/source/Solvers/solving_solar_cells.rst

@@ -97,7 +97,7 @@ The options available as well as the default values are:
         If the Boltzmann approximation should be used in the detailed balance solver. Any other choice will result in using the full Plank equation, which will be slower, in general.
 
 - :doc:`Depletion approximation solver options <depletion>`
-    - **da_mode** = 'bvp'
+    - **da_mode** = 'green'
         Selects the numerical approximation method for the drift-diffusion equation in the depletion approximation solver. Possible values are “bvp” for numerical solution using the `solve_bvp` method of the `scipy.integrate` module or 'green' for a semi-analytic solution using Green's functions. The latter is expected to be faster. 
 
 - :doc:`Poisson-drift diffusion solver options <DriftDiffusionUtilities>`

examples/MJ_solar_cell_external_reflectance.py

@@ -0,0 +1,125 @@
+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
+
+all_materials = []
+
+def this_dir_file(f):
+    return "data/" + f
+
+# 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)
+
+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)
+
+# BOTTOM CELL - Ge
+Ge = material("Ge")
+
+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
+
+# 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)
+solar_cell = SolarCell(
+    [
+        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'),
+        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'),
+        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,
+                                                  'da_mode': 'green'})
+
+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.legend()
+plt.ylim(0, 100)
+plt.ylabel('EQE (%)')
+plt.xlabel('Wavelength (nm)')
+plt.show()
+
+V = np.linspace(0, 3, 300)
+solar_cell_solver(solar_cell, 'iv', user_options={'voltages': V,
+                                                  'light_iv': True,
+                                                  'wavelength': wl,
+                                                  'da_mode': 'green'
+                                                  })
+
+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.legend()
+plt.ylim(0, 230)
+plt.xlim(0, 3)
+plt.ylabel('Current (A/m$^2$)')
+plt.xlabel('Voltage (V)')
+
+plt.show()