add another variant of implementing ohmic contact model

examples/Ex109_Traps.jl

@@ -159,6 +159,8 @@ function main(;n = 3, Plotter = PyPlot, plotting = false, verbose = false, test
     data.boundaryType[bregionAcceptor] = OhmicContact
     data.boundaryType[bregionDonor]    = OhmicContact
 
+    data.ohmicContactModel             = OhmicContactRobin
+
     ## Choose flux discretization scheme: ScharfetterGummel, ScharfetterGummelGraded,
     ## ExcessChemicalPotential, ExcessChemicalPotentialGraded, DiffusionEnhanced, GeneralizedSG
     data.fluxApproximation            .= ExcessChemicalPotential

src/ChargeTransport.jl

@@ -38,6 +38,8 @@ export OuterBoundaryModelType, OuterBoundaryModelType, InterfaceModelType
 export OhmicContact, SchottkyContact, SchottkyBarrierLowering, MixedOhmicSchottkyContact
 export InterfaceNone, InterfaceRecombination
 
+export OhmicContactModelType, OhmicContactDirichlet, OhmicContactRobin
+
 export ModelType, Transient, Stationary
 
 export FluxApproximationType

src/ct_datatypes.jl

@@ -81,6 +81,18 @@ Possible types of boundary models.
 """
 const BoundaryModelType  = Union{OuterBoundaryModelType, InterfaceModelType}
 
+##########################################################
+##########################################################
+
+abstract type OhmicContactDirichlet end
+
+abstract type OhmicContactRobin  end
+
+"""
+Possible mathematical types of ohmic contact boundary model.
+"""
+const OhmicContactModelType = Union{Type{OhmicContactDirichlet}, Type{OhmicContactRobin}}
+
 ##########################################################
 ##########################################################
 """

src/ct_physics.jl

@@ -299,6 +299,13 @@ breaction!(f, u, bnode, data) = breaction!(f, u, bnode, data, data.boundaryType[
 
 #################################################################################################
 
+breaction!(f, u, bnode, data, ::Type{OhmicContact}) = breaction!(f, u, bnode, data, data.calculationType)
+
+# in case of equilibrium conditions, we choose the initial ohmic contact boundary model
+breaction!(f, u, bnode, data, ::Type{InEquilibrium}) = breaction!(f, u, bnode, data, OhmicContactRobin)
+
+breaction!(f, u, bnode, data, ::Type{OutOfEquilibrium}) = breaction!(f, u, bnode, data, data.ohmicContactModel)
+
 """
 $(TYPEDSIGNATURES)
 
@@ -319,7 +326,7 @@ The boundary conditions for electrons and holes are dirichlet conditions, where
 
 with ``U`` as an applied voltage.
 """
-function breaction!(f, u, bnode, data, ::Type{OhmicContact})
+function breaction!(f, u, bnode, data, ::Type{OhmicContactRobin})
 
     params      = data.params
     paramsnodal = data.paramsnodal
@@ -389,6 +396,44 @@ function breaction!(f, u, bnode, data, ::Type{OhmicContact})
 
 end
 
+
+"""
+$(TYPEDSIGNATURES)
+
+Creates ohmic boundary conditions via Dirichlet BC for the electrostatic potential ``\\psi``
+
+``\\psi  = \\psi_0 + U``,
+
+where ``\\psi_0`` contains some given value and ``U`` is an applied voltage.
+
+``f[\\psi] =  -q/\\delta  \\sum_\\alpha{ z_\\alpha  (n_\\alpha - C_\\alpha) },``
+
+where ``C_\\alpha`` corresponds to some doping w.r.t. the species ``\\alpha``.
+
+The boundary conditions for electrons and holes are dirichlet conditions, where
+
+`` \\varphi_{\\alpha} = U.```
+
+"""
+function breaction!(f, u, bnode, data, ::Type{OhmicContactDirichlet})
+
+    params = data.params
+
+    # DA: we get here an issue with the allocation, if we pass into boundary_dirichlet! something which is not of type Int,
+    # as e.g. an AbstractQuantity
+    ipsi   = params.numberOfCarriers + 1 # data.index_psi
+    iphin  = data.bulkRecombination.iphin # integer index of φ_n
+    iphip  = data.bulkRecombination.iphip # integer index of φ_p
+
+    Δu     = params.contactVoltage[bnode.region] + data.contactVoltageFunction[bnode.region](bnode.time)
+    ψ0     = params.bψEQ[bnode.region]
+
+    boundary_dirichlet!(f, u, bnode, species = iphin, region=bnode.region, value=Δu)
+    boundary_dirichlet!(f, u, bnode, species = iphip, region=bnode.region, value=Δu)
+    boundary_dirichlet!(f, u, bnode, species = ipsi,  region=bnode.region, value=ψ0+Δu)
+
+end
+
 """
 $(TYPEDSIGNATURES)
 Creates Schottky boundary conditions. For the electrostatic potential we assume

src/ct_system.jl

@@ -358,10 +358,17 @@ mutable struct Params
     SchottkyBarrier              ::  Array{Float64,1}
 
     """
-    An array containing a contant value for the applied voltage.
+    An array containing a constant value for the applied voltage.
     """
     contactVoltage               ::  Array{Float64, 1}
 
+
+    """
+    An array containing a constant value for the electric potential
+    in case of Dirichlet boundary conditions.
+    """
+    bψEQ                         ::  Array{Float64, 1}
+
     ###############################################################
     ####                  number of carriers                   ####
     ###############################################################
@@ -715,6 +722,12 @@ mutable struct Data{TFuncs<:Function, TVoltageFunc<:Function, TGenerationData<:U
     """
     λ3                           ::  Float64
 
+    """
+    Possibility to change the implementation of the ohmic contact boundary model
+    for the electric potential (Dirichlet or Robin)
+    """
+    ohmicContactModel            :: OhmicContactModelType
+
     ###############################################################
     ####             Templates for DOS and BEE                 ####
     ###############################################################
@@ -832,6 +845,7 @@ function Params(grid, numberOfCarriers)
     ###############################################################
     params.SchottkyBarrier              = spzeros(Float64, numberOfBoundaryRegions)
     params.contactVoltage               = spzeros(Float64, numberOfBoundaryRegions)
+    params.bψEQ                         = spzeros(Float64, numberOfBoundaryRegions)
 
     ###############################################################
     ####                  number of carriers                   ####
@@ -855,6 +869,7 @@ function Params(grid, numberOfCarriers)
     params.bRecombinationSRHTrapDensity = spzeros(Float64, 2, numberOfBoundaryRegions)
     params.bRecombinationSRHLifetime    = spzeros(Float64, 2, numberOfBoundaryRegions)
     params.bDensityEQ                   = spzeros(Float64, 2, numberOfBoundaryRegions)
+
     ###############################################################
     ####        number of carriers x number of regions         ####
     ###############################################################
@@ -981,12 +996,13 @@ function Data(grid, numberOfCarriers; contactVoltageFunction = [zeroVoltage for
     ####                 Numerics information                  ####
     ###############################################################
     data.fluxApproximation                     = FluxApproximationType[ScharfetterGummel for i = 1:numberOfCarriers]
-    data.calculationType                       = OutOfEquilibrium     # do performances InEquilibrium or OutOfEquilibrium
-    data.modelType                             = Stationary        # indicates if we need additional time dependent part
-    data.generationModel                       = GenerationNone    # generation model
-    data.λ1                                    = 1.0               # λ1: embedding parameter for NLP
-    data.λ2                                    = 1.0               # λ2: embedding parameter for G
-    data.λ3                                    = 1.0               # λ3: embedding parameter for electro chemical reaction
+    data.calculationType                       = OutOfEquilibrium      # do performances InEquilibrium or OutOfEquilibrium
+    data.modelType                             = Stationary            # indicates if we need additional time dependent part
+    data.generationModel                       = GenerationNone        # generation model
+    data.λ1                                    = 1.0                   # λ1: embedding parameter for NLP
+    data.λ2                                    = 1.0                   # λ2: embedding parameter for G
+    data.λ3                                    = 1.0                   # λ3: embedding parameter for electro chemical reaction
+    data.ohmicContactModel                     = OhmicContactDirichlet # OhmicContactRobin also possible
 
     ###############################################################
     ####             Templates for DOS and BEE                 ####
@@ -1417,8 +1433,7 @@ gridplot(grid::ExtendableGrid; Plotter, kwargs...)                   = GridVisua
 """
 $(TYPEDSIGNATURES)
 
-Functions which sets for given charge carrier at a given boundary
-a given value.
+Functions which calculates the equilibrium solution in case of non-present fluxes and zero bias.
 
 """
 
@@ -1427,13 +1442,19 @@ function equilibrium_solve!(ctsys::System; control = VoronoiFVM.NewtonControl(),
     ctsys.fvmsys.physics.data.calculationType = InEquilibrium
     grid                                      = ctsys.fvmsys.grid
 
+    data        = ctsys.fvmsys.physics.data
+    params      = ctsys.fvmsys.physics.data.params
+    paramsnodal = ctsys.fvmsys.physics.data.paramsnodal
+    bnode       = grid[BFaceNodes]
+    ipsi        = data.index_psi
+
     # We set zero voltage for each charge carrier at all outer boundaries for equilibrium calculations.
     for ibreg ∈ grid[BFaceRegions]
         set_contact!(ctsys, ibreg, Δu = 0.0)
     end
 
     # initialize solution and starting vectors
-    if inival==nothing
+    if inival===nothing
         inival               = unknowns(ctsys)
         inival              .= 0.0
     end
@@ -1463,11 +1484,12 @@ function equilibrium_solve!(ctsys::System; control = VoronoiFVM.NewtonControl(),
 
     end
 
-    data        = ctsys.fvmsys.physics.data
-    params      = ctsys.fvmsys.physics.data.params
-    paramsnodal = ctsys.fvmsys.physics.data.paramsnodal
-    bnode       = grid[BFaceNodes]
-    ipsi        = data.index_psi
+    for ibreg ∈ grid[BFaceRegions]
+        # here we assume that in multidimensions, we receive a constant value of the electric potential at the boundary
+        # check for applications, where this is not the case
+        bψVal                   = view(sol[ipsi, :], subgrid(grid, [ibreg], boundary = true))[1]
+        params.bψEQ[ibreg] = bψVal
+    end
 
     # calculate equilibrium densities (especially needed for Schottky boundary conditions)
     for icc ∈ data.electricCarrierList