ViennaSHE 1.3.0
Free open-source semiconductor device simulator using spherical harmonics expansions techniques.
assemble_streaming.hpp
Go to the documentation of this file.
1#ifndef VIENNASHE_SHE_ASSEMBLE_STREAMING_HPP
2#define VIENNASHE_SHE_ASSEMBLE_STREAMING_HPP
3
4/* ============================================================================
5 Copyright (c) 2011-2022, Institute for Microelectronics,
6 Institute for Analysis and Scientific Computing,
7 TU Wien.
8
9 -----------------
10 ViennaSHE - The Vienna Spherical Harmonics Expansion Boltzmann Solver
11 -----------------
12
13 http://viennashe.sourceforge.net/
14
15 License: MIT (X11), see file LICENSE in the base directory
16=============================================================================== */
17
18
19// viennagrid
20#include "viennagrid/mesh/mesh.hpp"
21#include "viennagrid/algorithm/norm.hpp"
22#include "viennagrid/algorithm/volume.hpp"
23#include "viennagrid/algorithm/voronoi.hpp"
24
25// viennashe
28
32
33#include "viennashe/config.hpp"
38
42
43#include "viennashe/log/log.hpp"
44
49namespace viennashe
50{
51 namespace she
52 {
53
66 template <typename DeviceType,
67 typename SHEQuantity,
68 typename MatrixType,
69 typename VectorType,
70 typename CellType,
71 typename FacetType,
72 typename CouplingMatrixType>
74 viennashe::config const & conf,
75 SHEQuantity const & quan,
76 MatrixType & A,
77 VectorType & /*b*/,
78 CellType const & cell, FacetType const & facet, std::size_t index_H,
79 CouplingMatrixType const & coupling_matrix_diffusion,
80 CouplingMatrixType const & coupling_matrix_drift,
81 bool odd_assembly)
82 {
83 typedef typename viennagrid::result_of::point<CellType>::type PointType;
84
85 typename viennashe::config::dispersion_relation_type dispersion = conf.dispersion_relation(quan.get_carrier_type_id());
86
87 //const bool with_full_newton = (conf.nonlinear_solver().id() == viennashe::solvers::nonlinear_solver_ids::newton_nonlinear_solver);
88
89 CellType const *other_cell_ptr = util::get_other_cell_of_facet(device.mesh(), facet, cell);
90 if (!other_cell_ptr) return;
91
92 long row_index = odd_assembly ? quan.get_unknown_index(facet, index_H) : quan.get_unknown_index(cell, index_H);
93
94 if (row_index < 0) //no unknown here
95 return;
96
97 long col_index = odd_assembly ? quan.get_unknown_index(cell, index_H) : quan.get_unknown_index(facet, index_H);
98
99 if (col_index < 0) //other element does not carry an unknown, so nothing to do here
100 return;
101
102 PointType cell_connection = viennagrid::centroid(cell) - viennagrid::centroid(*other_cell_ptr);
103 const double connection_len = viennagrid::norm_2(cell_connection);
104 PointType cell_connection_normalized = cell_connection / viennagrid::norm(cell_connection);
105 PointType facet_unit_normal = viennashe::util::outer_cell_normal_at_facet(cell, facet);
106 const double weighted_interface_area = viennagrid::volume(facet) * std::fabs(viennagrid::inner_prod(facet_unit_normal, cell_connection_normalized));
107
108 long expansion_order_row = static_cast<long>(odd_assembly ? quan.get_expansion_order(facet, index_H) : quan.get_expansion_order(cell, index_H));
109 long expansion_order_column = static_cast<long>(odd_assembly ? quan.get_expansion_order(cell, index_H) : quan.get_expansion_order(facet, index_H));
110
111 double lower_kinetic_energy = lower_kinetic_energy_at_facet(device, quan, facet, index_H);
112 double upper_kinetic_energy = upper_kinetic_energy_at_facet(device, quan, facet, index_H);
113
118
119 //
120 // 'diffusion' term
121 //
122
123 double matrix_coeff_diffusion = 0;
124 switch (conf.she_discretization_type())
125 {
127 matrix_coeff_diffusion = integral_vZ(quan, dispersion, facet, index_H) * weighted_interface_area;
128 break;
130 matrix_coeff_diffusion = integral_v(quan, dispersion, facet, index_H) * weighted_interface_area;
131 break;
132 default: throw std::runtime_error("assemble_free_streaming_operator_on_box(): Unknown SHE discretization type!");
133 }
134
135 // write matrix contribution:
136 if (col_index < 0)
137 {
138 log::error() << "* assemble_free_streaming_operator(): Invalid column index detected while assembling diffusion operator. " << std::endl;
139 log::error<log_assemble_free_streaming_operator>() << " Diagnostics: " << std::endl << "index_H: " << index_H << std::endl;
140 log::error<log_assemble_free_streaming_operator>() << "cell: " << cell << std::endl;
141 log::error<log_assemble_free_streaming_operator>() << "row_index: " << row_index << std::endl;
142 log::error<log_assemble_free_streaming_operator>() << "facet: " << facet << std::endl;
143 log::error<log_assemble_free_streaming_operator>() << "odd_assembly: " << odd_assembly << std::endl;
144 log::error<log_assemble_free_streaming_operator>() << "DOF el2: " << col_index << std::endl;
145 log::error<log_assemble_free_streaming_operator>() << "other_cell: " << *other_cell_ptr << std::endl;
146 }
147
148 viennashe::util::add_block_matrix(A, std::size_t(row_index), std::size_t(col_index),
149 matrix_coeff_diffusion,
150 coupling_matrix_diffusion,
151 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
152 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col)
153 );
154
156 {
157 log::debug<log_assemble_free_streaming_operator>() << "* assemble_free_streaming_operator_on_box(): diffusion term contribution:" << std::endl;
158 log::debug<log_assemble_free_streaming_operator>() << " int_vZ * interface_area" << std::endl;
159 log::debug<log_assemble_free_streaming_operator>() << " at (" << row_index << ", " << col_index << ")" << std::endl;
160 log::debug<log_assemble_free_streaming_operator>() << " matrix_coeff_diffusion = " << matrix_coeff_diffusion << std::endl;
161 log::debug<log_assemble_free_streaming_operator>() << "weighted_interface_area = " << weighted_interface_area << std::endl;
162 log::debug<log_assemble_free_streaming_operator>() << " coupling_matrix_diffusion: " << coupling_matrix_diffusion << std::endl;
163 }
164
165 /*
166 if (with_full_newton)
167 {
168 double potential_offset = 0.05;
169 double energy_offset = potential_offset * viennashe::physics::constants::q;
170 double dint_vZ_dphi = (integral_vZ(lower_kinetic_energy + energy_offset,
171 upper_kinetic_energy + energy_offset,
172 dispersion)
173 - integral_vZ(lower_kinetic_energy,
174 upper_kinetic_energy,
175 dispersion)
176 ) / potential_offset;
177
178 double d_dphi = dint_vZ_dphi * interface_area;
179
180 if (odd_assembly)
181 {
182 // coupling with 'considered vertex' (edge-to-vertex coupling) is negative for electrons (energy, negative for holes
183 if (potential_index(vertex) >= 0)
184 viennashe::util::add_folded_block_matrix(A,
185 row_index, potential_index(vertex),
186 (assemble_electrons ? -0.5 : 0.5) * d_dphi,
187 coupling_matrix_diffusion,
188 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
189 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col),
190 current_guess, col_index);
191
192 // coupling with 'other' vertex (other vertex of edge) is negative for electrons, positive for holes:
193 if (potential_index(other_vertex) >= 0)
194 viennashe::util::add_folded_block_matrix(A,
195 row_index, potential_index(other_vertex),
196 (assemble_electrons ? 0.5 : -0.5) * d_dphi,
197 coupling_matrix_diffusion,
198 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
199 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col),
200 current_guess, col_index);
201 }
202 else
203 {
204 // coupling with 'this' vertex is negative for electrons, positive for holes
205 if (potential_index(vertex) >= 0)
206 viennashe::util::add_folded_block_matrix(A,
207 row_index, potential_index(vertex),
208 (assemble_electrons ? -0.5 : 0.5) * d_dphi,
209 coupling_matrix_diffusion,
210 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
211 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col),
212 current_guess, col_index);
213
214 // coupling with 'other' vertex is positive for electrons, negative for holes:
215 if (potential_index(other_vertex) >= 0)
216 viennashe::util::add_folded_block_matrix(A,
217 row_index, potential_index(other_vertex),
218 (assemble_electrons ? 0.5 : -0.5) * d_dphi,
219 coupling_matrix_diffusion,
220 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
221 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col),
222 current_guess, col_index);
223 }
224
225 // residual contribution:
226 viennashe::util::subtract_folded_block_vector(b,
227 row_index,
228 matrix_coeff_diffusion,
229 coupling_matrix_diffusion,
230 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
231 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col),
232 current_guess, col_index);
233 }
234 */
235
236 //
237 // 'drift' term
238 //
239
240 double Z = averaged_density_of_states(quan, dispersion, facet, index_H);
241 double int_Z_over_hk = integral_Z_over_hk(lower_kinetic_energy,
242 upper_kinetic_energy,
243 dispersion);
244 double force = force_on_facet_accessor()(quan, *other_cell_ptr, cell, index_H); //TODO: Here we need the global force, not just the projection on the edge
245 //force_along_edge(controller, other_vertex, vertex, index_H);
246
247 double matrix_coeff_drift = 0;
248 switch (conf.she_discretization_type())
249 {
251 matrix_coeff_drift = int_Z_over_hk * force * weighted_interface_area * connection_len / 2.0;
252 break;
254 matrix_coeff_drift = int_Z_over_hk / Z * force * weighted_interface_area * connection_len / 2.0;
255 break;
256 default: throw std::runtime_error("assemble_free_streaming_operator_on_box(): Unknown SHE discretization type!");
257 }
258
259
260 if (col_index < 0)
261 {
262 log::error() << "* assemble_free_streaming_operator(): Invalid column index detected while assembling drift operator. Diagnostics: " << std::endl;
263 log::error<log_assemble_free_streaming_operator>() << "index_H: " << index_H << std::endl;
264 log::error<log_assemble_free_streaming_operator>() << "cell: " << cell << std::endl;
265 log::error<log_assemble_free_streaming_operator>() << "row_index: " << row_index << std::endl;
266 log::error<log_assemble_free_streaming_operator>() << "facet: " << facet << std::endl;
267 log::error<log_assemble_free_streaming_operator>() << "odd_assembly: " << odd_assembly << std::endl;
268 log::error<log_assemble_free_streaming_operator>() << "DOF el2: " << col_index << std::endl;
269 log::error<log_assemble_free_streaming_operator>() << "other_cell: " << *other_cell_ptr << std::endl;
270 }
271
272 viennashe::util::add_block_matrix(A, std::size_t(row_index), std::size_t(col_index),
273 matrix_coeff_drift,
274 coupling_matrix_drift,
275 viennashe::math::spherical_harmonics_iterator(expansion_order_row, harmonics_it_id_row),
276 viennashe::math::spherical_harmonics_iterator(expansion_order_column, harmonics_it_id_col)
277 );
278
280 {
281 log::debug<log_assemble_free_streaming_operator>() << "* assemble_free_streaming_operator_on_box(): drift term contribution:" << std::endl;
282 log::debug<log_assemble_free_streaming_operator>() << " int_Z_over_hk / Z * force * interface_area * edge_len / 2.0" << std::endl;
283 log::debug<log_assemble_free_streaming_operator>() << " at (" << row_index << ", " << col_index << ")" << std::endl;
284 log::debug<log_assemble_free_streaming_operator>() << " int_Z_over_hk = " << int_Z_over_hk << std::endl;
285 log::debug<log_assemble_free_streaming_operator>() << " force = " << force << std::endl;
286 log::debug<log_assemble_free_streaming_operator>() << " weighted_interface_area = " << weighted_interface_area << std::endl;
287 log::debug<log_assemble_free_streaming_operator>() << " connection_len = " << connection_len << std::endl;
288 }
289
290
291 /*if (with_full_newton)
292 {
293 //We assume 'rate' to be independent of the potential, which is true for phonon scattering
294
295 double potential_offset = 0.05;
296 double energy_offset = potential_offset * viennashe::physics::constants::q;
297
298 double dint_Z_over_hk_dphi = (integral_Z_over_hk(lower_kinetic_energy + energy_offset, upper_kinetic_energy + energy_offset, dispersion)
299 - integral_Z_over_hk(lower_kinetic_energy, upper_kinetic_energy, dispersion)
300 ) / potential_offset;
301
302 double d_dphi = dint_Z_over_hk_dphi * force * interface_area * edge_len / 2.0; //Note: Ignoring dependence of force on the potential for the moment (term does not contribute for first-order SHE anyways)
303
304 std::vector<double> fold_vector(coupling_matrix_drift.size2());
305 if (assemble_electrons)
306 current_guess.fill(viennashe::ELECTRON_TYPE_ID, el2, controller.kinetic_energy(el2, index_H, viennashe::ELECTRON_TYPE_ID), index_H, fold_vector);
307 else
308 current_guess.fill(viennashe::HOLE_TYPE_ID, el2, controller.kinetic_energy(el2, index_H, viennashe::HOLE_TYPE_ID), index_H, fold_vector);
309
310 if (odd_assembly)
311 {
312 // coupling with 'considered vertex' (edge-to-vertex coupling) is positive for electrons, negative for holes
313 if (potential_index(vertex) >= 0)
314 viennashe::util::add_folded_block_matrix(A,
315 row_index, potential_index(vertex),
316 (assemble_electrons ? 0.5 : -0.5) * d_dphi,
317 coupling_matrix_drift,
318 HarmonicsIter1(expansion_order_row),
319 HarmonicsIter2(expansion_order_column),
320 fold_vector);
321
322 // coupling with 'other' vertex (other vertex of edge) is negative for electrons, positive for holes:
323 if (potential_index(other_vertex) >= 0)
324 viennashe::util::add_folded_block_matrix(A,
325 row_index, potential_index(other_vertex),
326 (assemble_electrons ? -0.5 : 0.5) * d_dphi,
327 coupling_matrix_drift,
328 HarmonicsIter1(expansion_order_row),
329 HarmonicsIter2(expansion_order_column),
330 fold_vector);
331 }
332 else
333 {
334 // coupling with 'this' vertex is negative for electrons, positive for holes
335 if (potential_index(vertex) >= 0)
336 viennashe::util::add_folded_block_matrix(A,
337 row_index, potential_index(vertex),
338 (assemble_electrons ? -1.0 : 1.0) * d_dphi,
339 coupling_matrix_drift,
340 HarmonicsIter1(expansion_order_row),
341 HarmonicsIter2(expansion_order_column),
342 fold_vector);
343
344 // coupling with 'other' vertex is positive for electrons, negative for holes:
345 if (potential_index(other_vertex) >= 0)
346 viennashe::util::add_folded_block_matrix(A,
347 row_index, potential_index(other_vertex),
348 (assemble_electrons ? 1.0 : -1.0) * d_dphi,
349 coupling_matrix_drift,
350 HarmonicsIter1(expansion_order_row),
351 HarmonicsIter2(expansion_order_column),
352 fold_vector);
353 }
354
355 // residual contribution:
356 viennashe::util::subtract_folded_block_vector(b,
357 row_index,
358 matrix_coeff,
359 coupling_matrix_drift,
360 HarmonicsIter1(expansion_order_row),
361 HarmonicsIter2(expansion_order_column),
362 fold_vector);
363 }*/
364
365 }
366
367
368
369 } //namespace she
370} //namespace viennashe
371
372#endif
Common routines for the assembly of SHE equations.
Writes a possibly scaled block matrix to the system matrix. Includes similar functionality for comput...
The main SHE configuration class. To be adjusted by the user for his/her needs.
Definition: config.hpp:124
viennashe::physics::dispersion_proxy dispersion_relation(viennashe::carrier_type_id ctype) const
Returns the dispersion relation for electrons.
Definition: config.hpp:266
she_discretization_type_id she_discretization_type() const
Definition: config.hpp:360
MeshT const & mesh() const
Returns the underlying mesh.
Definition: device.hpp:145
Defines the physical properties of a device, e.g. doping. This is the implementation for 2d and highe...
Definition: device.hpp:818
Iteration over all spherical harmonics indices up to order L, with increasing indices l and m.
A proxy object for a dispersion relation. Does NOT take ownership of the provided pointer!
Definition: dispersion.hpp:69
The SHE configuration class is defined here.
Contains the dispersion relations for different materials and different polarities.
Defines several filter functors for the device. A filter functor returns true if the supplied argumen...
Provides the SHE coupling matrices and computes higher-order coupling matrices if required....
Implementation of numerical integration routines.
A logging facility providing fine-grained control over logging in ViennaSHE.
Miscellaneous utilities.
logger< true > error()
Used to log errors. The logging level is logERROR.
Definition: log.hpp:301
harmonics_iteration_type
An enumeration of spherical harmonics types (all, even, odd) ViennaSHE supports.
Definition: forwards.h:159
@ EVEN_HARMONICS_ITERATION_ID
Definition: forwards.h:161
@ ODD_HARMONICS_ITERATION_ID
Definition: forwards.h:162
VectorType::value_type norm_2(VectorType const &v)
Definition: linalg_util.hpp:37
double integral_Z_over_hk(double energy_minus, double energy_plus, DispersionRelation const &dispersion)
Computes \int_a^b Z/(hbar * k) dE, where a and b are the energy boundaries, Z is the density of state...
double integral_v(SHEQuantity const &quan, viennashe::config::dispersion_relation_type const &dispersion, CellFacetType const &cell_facet, std::size_t index_H)
Computes \int v dH using a Simpson rule, where v is velocity.
double lower_kinetic_energy_at_facet(DeviceType const &device, SHEQuantity const &quan, FacetType const &facet, std::size_t index_H)
Returns the lower kinetic energy for the discretization box associated with the edge 'edge'.
double averaged_density_of_states(SHEQuantity const &quan, viennashe::config::dispersion_relation_type const &dispersion, CellFacetType const &cell_facet, std::size_t index_H)
Returns the density of states around a vertex or an edge at total energy specified by index_H....
double integral_vZ(SHEQuantity const &quan, viennashe::config::dispersion_relation_type const &dispersion, CellFacetType const &cell_facet, std::size_t index_H)
Computes \int v Z dH using a Simpson rule, where v is velocity and Z is the density of states.
double upper_kinetic_energy_at_facet(DeviceType const &device, SHEQuantity const &quan, FacetType const &facet, std::size_t index_H)
Returns the upper kinetic energy of the discretization box for the edge 'edge'.
void assemble_free_streaming_operator_on_box(DeviceType const &device, viennashe::config const &conf, SHEQuantity const &quan, MatrixType &A, VectorType &, CellType const &cell, FacetType const &facet, std::size_t index_H, CouplingMatrixType const &coupling_matrix_diffusion, CouplingMatrixType const &coupling_matrix_drift, bool odd_assembly)
Worker function for the assembly of the free streaming operator. Handles the assembly of both even an...
void add_block_matrix(SystemMatrixType &system_matrix, std::size_t row_index, std::size_t col_index, double prefactor, BlockMatrixType const &coupling_matrix, RowIndexIterator row_iter, ColumnIndexIterator const &col_iter_init, bool elastic_scattering_roundoff_error_prevention=false)
Writes a sub-matrix of a block matrix 'coupling_matrix' into the global system matrix,...
viennagrid::result_of::point< viennagrid::element< CellTagT, WrappedConfigT > >::type outer_cell_normal_at_facet(viennagrid::element< CellTagT, WrappedConfigT > const &cell, viennagrid::element< FacetTagT, WrappedConfigT > const &facet)
Returns the unit outer normal of a facet with respect to the provided cell.
CellT const * get_other_cell_of_facet(MeshT const &mesh, FacetT const &facet, CellT const &cell)
Helper function returning a const-pointer to the 'second cell' of a facet, or NULL if there is no sec...
Definition: misc.hpp:196
The main ViennaSHE namespace. All functionality resides inside this namespace.
Definition: accessors.hpp:40
@ SHE_DISCRETIZATION_EVEN_ODD_ORDER_GENERALIZED_DF
Definition: forwards.h:146
@ SHE_DISCRETIZATION_EVEN_ODD_ORDER_DF
Definition: forwards.h:145
Provides a number of fundamental constants. All constants in SI units.
Returns a few helper routines for computing physical quantities. To be replaced in the future.
Provides the exceptions used inside the viennashe::she namespace.
Implementation of spherical harmonics plus helper functions.
Defines the log keys used within the viennashe::she namespace.