?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
more fix on Ec/Ev.
[gss-tcad.git] / src / solver / mix2 / semiequ2mix.cc
blob4f4d597a4bdf90e90ab0806bf0de8fba7707d4ea
1 /*****************************************************************************/
2 /* 8888888 88888888 88888888 */
3 /* 8 8 8 */
4 /* 8 8 8 */
5 /* 8 88888888 88888888 */
6 /* 8 8888 8 8 */
7 /* 8 8 8 8 */
8 /* 888888 888888888 888888888 */
9 /* */
10 /* A Two-Dimensional General Purpose Semiconductor Simulator. */
11 /* */
12 /* GSS 0.4x */
13 /* Last update: May 11, 2007 */
14 /* */
15 /* Gong Ding gdiso@ustc.edu */
16 /* Xuan Chun xiaomoyu505@163.com */
17 /* NINT, No.69 P.O.Box, Xi'an City, China */
18 /* */
19 /*****************************************************************************/
20
21 #include "zonedata.h"
22 #include "vec4.h"
23
24
25 //-----------------------------------------------------------------------------
26 // special process of electrode bcs for mix type simulation
27 //-----------------------------------------------------------------------------
28
29
30 void SMCZone::F2E_mix_ddm_ombc_segment(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
31 {
32 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
33 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
34 int z = zone_index;
35 int size = pzone->davcell.size();
36 PetscScalar e = mt->e;
37 PetscScalar kb = mt->kb;
38 PetscScalar Na = aux[i].Na;
39 PetscScalar Nd = aux[i].Nd;
40 PetscScalar Vi = x[zofs[z]+4*i+0]; //potential of node i
41 PetscScalar ni = x[zofs[z]+4*i+1]; //electron density of node i
42 PetscScalar pi = x[zofs[z]+4*i+2]; //hole density of node i
43 PetscScalar Ti = fabs(x[zofs[z]+4*i+3]);
44 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
45 PetscScalar nie = mt->band->nie(fs[i].T);
46 PetscScalar Nc = mt->band->Nc(fs[i].T);
47 PetscScalar Nv = mt->band->Nv(fs[i].T);
48 PetscScalar H = 0;
49 PetscScalar electron_density,hole_density;
50
51 if(Fermi) //Fermi
52 {
53 PetscScalar Ec = -(e*Vi + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) + kb*fs[i].T*log(aux[i].Nc));
54 PetscScalar Ev = -(e*Vi + aux[i].affinity - mt->band->EgNarrowToEv(fs[i].T) - kb*fs[i].T*log(aux[i].Nv)+ aux[i].Eg);
55 PetscScalar phin = pbc->Vapp;
56 PetscScalar phip = pbc->Vapp;
57 PetscScalar etan = (-e*phin-Ec)/kb/fs[i].T;
58 PetscScalar etap = (Ev+e*phip)/kb/fs[i].T;
59 f[zofs[z]+4*i+0] = Nc*fermi_half(etan) - Nv*fermi_half(etap)+ Na - Nd;
60 f[zofs[z]+4*i+1] = ni - Nc*fermi_half(etan);
61 f[zofs[z]+4*i+2] = pi - Nv*fermi_half(etap);
62 }
63 else
64 {
65 f[zofs[z]+4*i+0] = Vi - kb*fs[i].T/mt->e*asinh((Nd-Na)/(2*nie)) + mt->kb*fs[i].T/mt->e*log(Nc/nie)
66 + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) - pbc->Vapp;
67
68 if(Na>Nd) //p-type
69 {
70 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
71 electron_density = nie*nie/hole_density;
72 }
73 else //n-type
74 {
75 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
76 hole_density = nie*nie/electron_density;
77 }
78 f[zofs[z]+4*i+1] = x[zofs[z]+4*i+1] - electron_density; //electron density
79 f[zofs[z]+4*i+2] = x[zofs[z]+4*i+2] - hole_density; //hole density
80 }
81 PetscScalar grad_T=0;
82 for(int j=0;j<pcell->nb_num;j++)
83 {
84 int nb = pcell->nb_array[j];
85 PetscScalar Tj = x[zofs[zone_index]+4*nb+3]; //lattice temperature of nb node
86 PetscScalar T_mid = 0.5*(Tj+Ti);
87 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
88 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
89 grad_T += kapa*pcell->elen[j]*dTdx_mid/pcell->area;
90 if(j==0||j==pcell->nb_num-1)
91 {
92 PetscScalar h = pbc->heat_transfer;
93 PetscScalar r = h*pbc->T_external;
94 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj))/pcell->area;
95 }
96 }
97 f[zofs[zone_index]+4*i+3] = grad_T + H;
98
99 if(ODE_F.TimeDependent==true)
100 {
101 //second order
102 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
103 {
104 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
105 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
106 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
107 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
108 }
109 else //first order
110 {
111 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
112 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
113 }
114 }
115
116 }
117
118
119 void SMCZone::F2E_mix_ddm_ombc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
120 ElZone *pz, int n)
121 {
122 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
123 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
124 int z = zone_index;
125 int equ_num = 4;
126 int size = pzone->davcell.size();
127 PetscScalar e = mt->e;
128 PetscScalar kb = mt->kb;
129 PetscScalar Na = aux[i].Na;
130 PetscScalar Nd = aux[i].Nd;
131 PetscScalar Vi = x[zofs[z]+4*i+0]; //potential of node i
132 PetscScalar ni = x[zofs[z]+4*i+1]; //electron density of node i
133 PetscScalar pi = x[zofs[z]+4*i+2]; //hole density of node i
134 PetscScalar Ti = fabs(x[zofs[z]+4*i+3]);
135 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
136 PetscScalar nie = mt->band->nie(fs[i].T);
137 PetscScalar Nc = mt->band->Nc(fs[i].T);
138 PetscScalar Nv = mt->band->Nv(fs[i].T);
139 PetscScalar H = 0;
140 PetscScalar electron_density,hole_density;
141
142 if(Fermi) //Fermi
143 {
144 PetscScalar Ec = -(e*Vi + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) + kb*fs[i].T*log(aux[i].Nc));
145 PetscScalar Ev = -(e*Vi + aux[i].affinity - mt->band->EgNarrowToEv(fs[i].T) - kb*fs[i].T*log(aux[i].Nv)+ aux[i].Eg);
146 PetscScalar phin = pbc->Vapp;
147 PetscScalar phip = pbc->Vapp;
148 PetscScalar etan = (-e*phin-Ec)/kb/fs[i].T;
149 PetscScalar etap = (Ev+e*phip)/kb/fs[i].T;
150 f[zofs[z]+4*i+0] = Nc*fermi_half(etan) - Nv*fermi_half(etap)+ Na - Nd;
151 f[zofs[z]+4*i+1] = ni - Nc*fermi_half(etan);
152 f[zofs[z]+4*i+2] = pi - Nv*fermi_half(etap);
153 }
154 {
155 f[zofs[z]+4*i+0] = Vi - kb*fs[i].T/mt->e*asinh((Nd-Na)/(2*nie)) + mt->kb*fs[i].T/mt->e*log(Nc/nie)
156 + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) - pbc->Vapp;
157
158 if(Na>Nd) //p-type
159 {
160 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
161 electron_density = nie*nie/hole_density;
162 }
163 else //n-type
164 {
165 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
166 hole_density = nie*nie/electron_density;
167 }
168 f[zofs[z]+4*i+1] = x[zofs[z]+4*i+1] - electron_density; //electron density
169 f[zofs[z]+4*i+2] = x[zofs[z]+4*i+2] - hole_density; //hole density
170 }
171 PetscScalar grad_T=0;
172 for(int j=0;j<pcell->nb_num;j++)
173 {
174 int nb = pcell->nb_array[j];
175 PetscScalar Tj = x[zofs[zone_index]+4*nb+3]; //lattice temperature of nb node
176 PetscScalar T_mid = 0.5*(Tj+Ti);
177 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
178 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
179 grad_T += kapa*pcell->elen[j]*dTdx_mid;
180 }
181 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
182 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
183 for(int j=0;j<ncell->nb_num;j++)
184 {
185 int nb = ncell->nb_array[j];
186 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
187 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
188 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
189 }
190 f[zofs[zone_index]+4*i+3] = grad_T + H;
191
192 if(ODE_F.TimeDependent==true)
193 {
194 //second order
195 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
196 {
197 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
198 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
199 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
200 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
201 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
202 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
203 }
204 else //first order
205 {
206 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
207 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
208 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
209 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
210 }
211 }
212
213 }
214
215
216 void SMCZone::F2E_mix_ddm_stkbc_segment(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc)
217 {
218 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
219 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
220 int z = zone_index;
221 int size = pzone->davcell.size();
222 PetscScalar kb = mt->kb;
223 PetscScalar e = mt->e;
224 PetscScalar Vi = x[zofs[zone_index]+4*i+0]; //potential of node i
225 PetscScalar ni = x[zofs[zone_index]+4*i+1]; //electron density of node i
226 PetscScalar pi = x[zofs[zone_index]+4*i+2]; //hole density of node i
227 PetscScalar Ti = x[zofs[zone_index]+4*i+3]; //lattice temperature of node i
228 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
229
230 //Schotty Barrier Lowerring
231 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
232 //potential
233 f[zofs[z]+4*i+0] = x[zofs[z]+4*i+0] + pbc->WorkFunction - deltaVB - pbc->Vapp;
234
235 PetscScalar Fn=0,Fp=0,grad_T=0;
236 for(int j=0;j<pcell->nb_num;j++)
237 {
238 int nb = pcell->nb_array[j];
239 PetscScalar Tj = x[zofs[zone_index]+4*nb+3]; //lattice temperature of nb node
240 if(j==0||j==pcell->nb_num-1)
241 {
242 Fn += mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[i].affinity - deltaVB)*0.5*pcell->ilen[j];
243 Fp += mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[i].affinity + deltaVB)*0.5*pcell->ilen[j];
244 PetscScalar h = pbc->heat_transfer;
245 PetscScalar r = h*pbc->T_external;
246 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj));
247 }
248 }
249
250 f[zofs[zone_index]+4*i+1] = (f[zofs[zone_index]+4*i+1]+Fn)/pcell->area;
251 f[zofs[zone_index]+4*i+2] = (f[zofs[zone_index]+4*i+2]-Fp)/pcell->area;
252 f[zofs[zone_index]+4*i+3] = (f[zofs[zone_index]+4*i+3]+grad_T)/pcell->area;
253
254 if(ODE_F.TimeDependent==true)
255 {
256 //second order
257 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
258 {
259 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
260 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
261 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
262 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
263 f[zofs[zone_index]+4*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
264 f[zofs[zone_index]+4*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
265 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
266 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
267 }
268 else //first order
269 {
270 f[zofs[zone_index]+4*i+1] += -(ni-fs[i].n)/ODE_F.dt;
271 f[zofs[zone_index]+4*i+2] += -(pi-fs[i].p)/ODE_F.dt;
272 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
273 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
274 }
275 }
276 }
277
278
279 void SMCZone::F2E_mix_ddm_stkbc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc,
280 ElZone *pz, int n)
281 {
282 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
283 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
284 int z = zone_index;
285 int equ_num = 4;
286 int size = pzone->davcell.size();
287 PetscScalar kb = mt->kb;
288 PetscScalar e = mt->e;
289 PetscScalar Vi = x[zofs[zone_index]+4*i+0]; //potential of node i
290 PetscScalar ni = x[zofs[zone_index]+4*i+1]; //electron density of node i
291 PetscScalar pi = x[zofs[zone_index]+4*i+2]; //hole density of node i
292 PetscScalar Ti = x[zofs[zone_index]+4*i+3]; //lattice temperature of node i
293 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
294
295 //Schotty Barrier Lowerring
296 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
297 //potential
298 f[zofs[z]+4*i+0] = x[zofs[z]+4*i+0] + pbc->WorkFunction - deltaVB - pbc->Vapp;
299
300 PetscScalar Fn=0,Fp=0,grad_T=0;
301 for(int j=0;j<pcell->nb_num;j++)
302 {
303 int nb = pcell->nb_array[j];
304 if(j==0||j==pcell->nb_num-1)
305 {
306 Fn += mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[i].affinity - deltaVB)*0.5*pcell->ilen[j];
307 Fp += mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[i].affinity + deltaVB)*0.5*pcell->ilen[j];
308 }
309 }
310 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
311 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
312 for(int j=0;j<ncell->nb_num;j++)
313 {
314 int nb = ncell->nb_array[j];
315 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
316 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
317 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
318 }
319 f[zofs[zone_index]+4*i+1] = (f[zofs[zone_index]+4*i+1]+Fn)/pcell->area;
320 f[zofs[zone_index]+4*i+2] = (f[zofs[zone_index]+4*i+2]-Fp)/pcell->area;
321 f[zofs[zone_index]+4*i+3] = (f[zofs[zone_index]+4*i+3]+grad_T);
322
323 if(ODE_F.TimeDependent==true)
324 {
325 //second order
326 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
327 {
328 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
329 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
330 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
331 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
332 f[zofs[zone_index]+4*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
333 f[zofs[zone_index]+4*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
334 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
335 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
336 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
337 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
338 }
339 else //first order
340 {
341 f[zofs[zone_index]+4*i+1] += -(ni-fs[i].n)/ODE_F.dt;
342 f[zofs[zone_index]+4*i+2] += -(pi-fs[i].p)/ODE_F.dt;
343 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
344 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
345 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
346 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
347 }
348 }
349 }
350
351
352 void SMCZone::F2E_mix_ddm_insulator_gate(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
353 {
354 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
355 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
356
357 int size = pzone->davcell.size();
358 PetscScalar kb = mt->kb;
359 PetscScalar e = mt->e;
360 PetscScalar Vi = x[zofs[zone_index]+4*i+0]; //potential of node i
361 PetscScalar ni = x[zofs[zone_index]+4*i+1]; //electron density of node i
362 PetscScalar pi = x[zofs[zone_index]+4*i+2]; //hole density of node i
363 PetscScalar Ti = x[zofs[zone_index]+4*i+3]; //lattice temperature of node i
364 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
365 PetscScalar L = 0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
366 PetscScalar grad_P = 0, grad_T = 0;
367
368 for(int j=0;j<pcell->nb_num;j++)
369 {
370 int nb = pcell->nb_array[j];
371 PetscScalar Vj = x[zofs[zone_index]+4*nb+0]; //potential of nb node
372 PetscScalar Tj = x[zofs[zone_index]+4*nb+3]; //lattice temperature of nb node
373
374 if(j==0||j==pcell->nb_num-1)
375 {
376 PetscScalar vgate = pbc->Vapp - pbc->WorkFunction;
377 PetscScalar q = e*pbc->QF; //sigma is the surface change density
378 PetscScalar Thick = pbc->Thick;
379 PetscScalar eps_ox = mt->eps0*pbc->eps;
380 PetscScalar eps = aux[i].eps;
381 PetscScalar r=q/eps + eps_ox/eps/Thick*vgate;
382 PetscScalar s=eps_ox/eps/Thick;
383 grad_P += 0.5*pcell->ilen[j]*(r-0.25*s*(3*Vi+Vj));
384
385 PetscScalar h = pbc->heat_transfer;
386 grad_T += 0.5*pcell->ilen[j]*(h*pbc->T_external-0.25*h*(3*Ti+Tj));
387 }
388 }
389
390 f[zofs[zone_index]+4*i+0] = (f[zofs[zone_index]+4*i+0]+grad_P)/pcell->area
391 + e/aux[i].eps*((pi-aux[i].Na)+(aux[i].Nd-ni));
392 f[zofs[zone_index]+4*i+1] = f[zofs[zone_index]+4*i+1]/pcell->area;
393 f[zofs[zone_index]+4*i+2] = f[zofs[zone_index]+4*i+2]/pcell->area;
394 f[zofs[zone_index]+4*i+3] = (f[zofs[zone_index]+4*i+3]+grad_T)/pcell->area;
395
396
397 if(ODE_F.TimeDependent==true)
398 {
399 //second order
400 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
401 {
402 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
403 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
404 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
405 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
406 f[zofs[zone_index]+4*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
407 f[zofs[zone_index]+4*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
408 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
409 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
410 }
411 else //first order
412 {
413 f[zofs[zone_index]+4*i+1] += -(ni-fs[i].n)/ODE_F.dt;
414 f[zofs[zone_index]+4*i+2] += -(pi-fs[i].p)/ODE_F.dt;
415 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
416 f[zofs[zone_index]+4*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
417 }
418 }
419 }
420
421
422 //------------------------------------------------------------------------------------------
423
424 void SMCZone::J2E_mix_ddm_ombc_segment(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
425 {
426 Mat4 A;
427 PetscInt index[4],col[4];
428 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
429 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
430 int z = zone_index;
431 PetscScalar kb = mt->kb;
432 PetscScalar e = mt->e;
433 PetscScalar Ti = x[zofs[z]+4*i+3]; //lattice temperature of node i
434 PetscScalar Vt = kb*Ti/e;
435 PetscScalar area = pcell->area;
436 PetscScalar d_grad_T_dTi = 0;
437 index[0] = zofs[z]+4*i+0;
438 index[1] = zofs[z]+4*i+1;
439 index[2] = zofs[z]+4*i+2;
440 index[3] = zofs[z]+4*i+3;
441 for(int j=0;j<pcell->nb_num;j++)
442 {
443 int nb = pcell->nb_array[j];
444 PetscScalar Tj = x[zofs[z]+4*nb+3]; //lattice temperature of nb node
445 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
446 //-------------------------------------
447 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j]/area;
448 //-------------------------------------
449 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]/area; //dfun(3)/dT(r)
450 if(j==0||j==pcell->nb_num-1)
451 {
452 PetscScalar h = pbc->heat_transfer;
453 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/area;
454 d_grad_T_dTj += -0.5*pcell->ilen[j]*0.25*h/area;
455 }
456 MatSetValue(*jac,zofs[z]+4*i+3,zofs[z]+4*nb+3,d_grad_T_dTj,INSERT_VALUES);
457 }
458 if(Fermi) //Fermi
459 {
460 PetscScalar Vt = kb*fs[i].T/e;
461 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
462 PetscScalar Nc = mt->band->Nc(fs[i].T);
463 PetscScalar Nv = mt->band->Nv(fs[i].T);
464 PetscScalar Vi = x[zofs[zone_index]+4*i+0]; //potential of node i
465 PetscScalar Ec = -(e*Vi + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) + kb*fs[i].T*log(aux[i].Nc));
466 PetscScalar Ev = -(e*Vi + aux[i].affinity - mt->band->EgNarrowToEv(fs[i].T) - kb*fs[i].T*log(aux[i].Nv)+ aux[i].Eg);
467 PetscScalar phin = pbc->Vapp;
468 PetscScalar phip = pbc->Vapp;
469 PetscScalar etan = (-e*phin-Ec)/kb/fs[i].T;
470 PetscScalar etap = (Ev+e*phip)/kb/fs[i].T;
471 PetscScalar detan_dVi = 1.0/Vt;
472 PetscScalar detap_dVi = -1.0/Vt;
473 Set_Mat4_zero(A);
474 A.m[0] = Nc*fermi_mhalf(etan)*detan_dVi - Nv*fermi_mhalf(etap)*detap_dVi;
475 A.m[4] = -Nc*fermi_mhalf(etan)*detan_dVi;
476 A.m[5] = 1;
477 A.m[8] = -Nv*fermi_mhalf(etap)*detap_dVi;
478 A.m[10] = 1;
479 A.m[15] = d_grad_T_dTi;
480 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
481 }
482 else //Boltzmann
483 {
484 Set_Mat4_I(A);
485 A.m[15] = d_grad_T_dTi; //dfun(3)/dT(i)
486 }
487 if(ODE_F.TimeDependent==true)
488 {
489 //second order
490 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
491 {
492 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
493 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
494 A.m[15] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
495 }
496 else //first order
497 {
498 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
499 A.m[15] += -aux[i].density*HeatCapacity/ODE_F.dt;
500 }
501 }
502
503 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
504 }
505
506
507 void SMCZone::J2E_mix_ddm_ombc_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc,
508 ElZone *pz, int n)
509 {
510 Mat4 A;
511 PetscInt index[4],col[4];
512 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
513 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
514 int z = zone_index;
515 PetscScalar kb = mt->kb;
516 PetscScalar e = mt->e;
517 PetscScalar Ti = x[zofs[z]+4*i+3]; //lattice temperature of node i
518 PetscScalar Vt = kb*Ti/e;
519 PetscScalar area = pcell->area;
520 PetscScalar d_grad_T_dTi = 0;
521 index[0] = zofs[z]+4*i+0;
522 index[1] = zofs[z]+4*i+1;
523 index[2] = zofs[z]+4*i+2;
524 index[3] = zofs[z]+4*i+3;
525 for(int j=0;j<pcell->nb_num;j++)
526 {
527 int nb = pcell->nb_array[j];
528 PetscScalar Tj = x[zofs[z]+4*nb+3]; //lattice temperature of nb node
529 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
530 //-------------------------------------
531 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j];
532 //-------------------------------------
533 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]; //dfun(3)/dT(r)
534 MatSetValue(*jac,zofs[z]+4*i+3,zofs[z]+4*nb+3,d_grad_T_dTj,INSERT_VALUES);
535 }
536 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
537 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
538 for(int j=0;j<ncell->nb_num;j++)
539 {
540 int nb = ncell->nb_array[j];
541 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
542 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
543 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
544 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
545 MatSetValue(*jac,zofs[z]+4*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
546 }
547
548 if(Fermi) //Fermi
549 {
550 PetscScalar Vt = kb*fs[i].T/e;
551 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
552 PetscScalar Nc = mt->band->Nc(fs[i].T);
553 PetscScalar Nv = mt->band->Nv(fs[i].T);
554 PetscScalar Vi = x[zofs[zone_index]+4*i+0]; //potential of node i
555 PetscScalar Ec = -(e*Vi + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) + kb*fs[i].T*log(aux[i].Nc));
556 PetscScalar Ev = -(e*Vi + aux[i].affinity - mt->band->EgNarrowToEv(fs[i].T) - kb*fs[i].T*log(aux[i].Nv)+ aux[i].Eg);
557 PetscScalar phin = pbc->Vapp;
558 PetscScalar phip = pbc->Vapp;
559 PetscScalar etan = (-e*phin-Ec)/kb/fs[i].T;
560 PetscScalar etap = (Ev+e*phip)/kb/fs[i].T;
561 PetscScalar detan_dVi = 1.0/Vt;
562 PetscScalar detap_dVi = -1.0/Vt;
563 Set_Mat4_zero(A);
564 A.m[0] = Nc*fermi_mhalf(etan)*detan_dVi - Nv*fermi_mhalf(etap)*detap_dVi;
565 A.m[4] = -Nc*fermi_mhalf(etan)*detan_dVi;
566 A.m[5] = 1;
567 A.m[8] = -Nv*fermi_mhalf(etap)*detap_dVi;
568 A.m[10] = 1;
569 A.m[15] = d_grad_T_dTi;
570 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
571 }
572 else //Boltzmann
573 {
574 Set_Mat4_I(A);
575 A.m[15] = d_grad_T_dTi; //dfun(3)/dT(i)
576 }
577
578 if(ODE_F.TimeDependent==true)
579 {
580 //second order
581 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
582 {
583 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
584 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
585 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
586 A.m[15] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
587 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
588 }
589 else //first order
590 {
591 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
592 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
593 A.m[15] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
594 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
595 }
596 }
597
598 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
599 }
600
601
602 void SMCZone::J2E_mix_ddm_stkbc_segment(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
603 {
604 Mat4 A;
605 PetscInt index[4],col[4];
606 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
607 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
608 int z = zone_index;
609 int size = pzone->davcell.size();
610
611 PetscScalar kb = mt->kb;
612 PetscScalar e = mt->e;
613 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
614 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
615 PetscScalar ni = x[zofs[z]+4*i+1]; //electron density of node i
616 PetscScalar pi = x[zofs[z]+4*i+2]; //hole density of node i
617 PetscScalar Ti = x[zofs[z]+4*i+3]; //lattice temperature of node i
618
619 PetscScalar area = pcell->area;
620 PetscScalar d_grad_T_dTi = 0;
621 PetscScalar d_Fn_dni = 0;
622 PetscScalar d_Fn_dTi = 0;
623 PetscScalar d_Fp_dpi = 0;
624 PetscScalar d_Fp_dTi = 0;
625 //--------------------------------
626 index[0] = zofs[z]+4*i+0;
627 index[1] = zofs[z]+4*i+1;
628 index[2] = zofs[z]+4*i+2;
629 index[3] = zofs[z]+4*i+3;
630 //--------------------------------
631 for(int j=0;j<pcell->nb_num;j++)
632 {
633 int nb = pcell->nb_array[j];
634 PetscScalar nj = x[zofs[z]+4*nb+1]; //electron density of nb node
635 PetscScalar pj = x[zofs[z]+4*nb+2]; //hole density of nb node
636 PetscScalar Tj = x[zofs[z]+4*nb+3]; //lattice temperature of nb node
637 col[0] = zofs[z]+4*nb+0;
638 col[1] = zofs[z]+4*nb+1;
639 col[2] = zofs[z]+4*nb+2;
640 col[3] = zofs[z]+4*nb+3;
641 MatGetValues(*jtmp,4,index,4,col,A.m);
642 A.m[0] = 0.0;
643 A=A/area;
644 if(j==0||j==pcell->nb_num-1)
645 {
646 d_Fn_dni += mt->band->pdSchottyJsn_pdn (ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
647 d_Fn_dTi += mt->band->pdSchottyJsn_pdTl(ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
648 d_Fp_dpi += mt->band->pdSchottyJsp_pdp (pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
649 d_Fp_dTi += mt->band->pdSchottyJsp_pdTl(pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
650 PetscScalar h = pbc->heat_transfer;
651 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
652 A.m[15] += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
653 }
654 MatSetValues(*jac,4,index,4,col,A.m,INSERT_VALUES);
655 }
656 //fun(0) is the poisson's equation
657 //fun(1) is the continuous equation of electron
658 //fun(2) is the continuous equation of hole
659 MatGetValues(*jtmp,4,index,4,index,A.m);
660 A=A/area;
661
662 A.m[0] = 1;
663
664 A.m[5] += d_Fn_dni; //dfun(1)/dn(i)
665 A.m[7] += d_Fn_dTi;
666
667 A.m[10]+= -d_Fp_dpi; //dfun(2)/dp(i)
668 A.m[11]+= -d_Fp_dTi;
669
670 A.m[15]+= d_grad_T_dTi; //dfun(3)/dT(i)
671
672 if(ODE_F.TimeDependent==true)
673 {
674 //second order
675 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
676 {
677 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
678 A.m[5] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
679 A.m[10] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
680 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
681 A.m[15] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
682 }
683 else //first order
684 {
685 A.m[5] += -1/ODE_F.dt;
686 A.m[10] += -1/ODE_F.dt;
687 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
688 A.m[15] += -aux[i].density*HeatCapacity/ODE_F.dt;
689 }
690 }
691 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
692 }
693
694
695 void SMCZone::J2E_mix_ddm_stkbc_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc,
696 ElZone *pz, int n)
697 {
698 Mat4 A;
699 PetscInt index[4],col[4];
700 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
701 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
702 int z = zone_index;
703 int equ_num = 4;
704 int size = pzone->davcell.size();
705 PetscScalar area = pcell->area;
706 Vec4 scale;
707 Set_Vec4(scale,area,area,area,1.0);
708 PetscScalar kb = mt->kb;
709 PetscScalar e = mt->e;
710 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
711 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
712 PetscScalar ni = x[zofs[z]+4*i+1]; //electron density of node i
713 PetscScalar pi = x[zofs[z]+4*i+2]; //hole density of node i
714 PetscScalar Ti = x[zofs[z]+4*i+3]; //lattice temperature of node i
715
716 PetscScalar d_grad_T_dTi = 0;
717 PetscScalar d_Fn_dni = 0;
718 PetscScalar d_Fn_dTi = 0;
719 PetscScalar d_Fp_dpi = 0;
720 PetscScalar d_Fp_dTi = 0;
721 //--------------------------------
722 index[0] = zofs[z]+4*i+0;
723 index[1] = zofs[z]+4*i+1;
724 index[2] = zofs[z]+4*i+2;
725 index[3] = zofs[z]+4*i+3;
726 //--------------------------------
727 for(int j=0;j<pcell->nb_num;j++)
728 {
729 int nb = pcell->nb_array[j];
730 PetscScalar nj = x[zofs[z]+4*nb+1]; //electron density of nb node
731 PetscScalar pj = x[zofs[z]+4*nb+2]; //hole density of nb node
732 PetscScalar Tj = x[zofs[z]+4*nb+3]; //lattice temperature of nb node
733 col[0] = zofs[z]+4*nb+0;
734 col[1] = zofs[z]+4*nb+1;
735 col[2] = zofs[z]+4*nb+2;
736 col[3] = zofs[z]+4*nb+3;
737 MatGetValues(*jtmp,4,index,4,col,A.m);
738 A.m[0] = 0.0;
739 A=A/scale;
740 if(j==0||j==pcell->nb_num-1)
741 {
742 d_Fn_dni += mt->band->pdSchottyJsn_pdn (ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
743 d_Fn_dTi += mt->band->pdSchottyJsn_pdTl(ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
744 d_Fp_dpi += mt->band->pdSchottyJsp_pdp (pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
745 d_Fp_dTi += mt->band->pdSchottyJsp_pdTl(pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
746 }
747 MatSetValues(*jac,4,index,4,col,A.m,INSERT_VALUES);
748 }
749 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
750 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
751 for(int j=0;j<ncell->nb_num;j++)
752 {
753 int nb = ncell->nb_array[j];
754 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
755 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
756 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
757 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
758 MatSetValue(*jac,zofs[z]+4*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
759 }
760 //fun(0) is the poisson's equation
761 //fun(1) is the continuous equation of electron
762 //fun(2) is the continuous equation of hole
763 MatGetValues(*jtmp,4,index,4,index,A.m);
764 A=A/scale;
765
766 A.m[0] = 1;
767
768 A.m[5] += d_Fn_dni; //dfun(1)/dn(i)
769 A.m[7] += d_Fn_dTi;
770
771 A.m[10]+= -d_Fp_dpi; //dfun(2)/dp(i)
772 A.m[11]+= -d_Fp_dTi;
773
774 A.m[15]+= d_grad_T_dTi; //dfun(3)/dT(i)
775
776 if(ODE_F.TimeDependent==true)
777 {
778 //second order
779 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
780 {
781 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
782 A.m[5] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
783 A.m[10] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
784 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
785 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
786 A.m[15] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
787 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
788 }
789 else //first order
790 {
791 A.m[5] += -1/ODE_F.dt;
792 A.m[10] += -1/ODE_F.dt;
793 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
794 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
795 A.m[15] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
796 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
797 }
798 }
799 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
800 }
801
802
803 void SMCZone::J2E_mix_ddm_insulator_gate(int i,PetscScalar *x, Mat *jac, Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
804 {
805 Mat4 A;
806 PetscInt index[4],col[4];
807 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
808 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
809 int z = zone_index;
810 int size = pzone->davcell.size();
811 PetscScalar kb = mt->kb;
812 PetscScalar e = mt->e;
813 PetscScalar Vi = x[zofs[z]+4*i+0]; //potential of node i
814 PetscScalar ni = x[zofs[z]+4*i+1]; //electron density of node i
815 PetscScalar pi = x[zofs[z]+4*i+2]; //hole density of node i
816 PetscScalar Ti = x[zofs[z]+4*i+3]; //lattice temperature of node i
817
818 PetscScalar area = pcell->area;
819 PetscScalar L = 0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
820 PetscScalar d_grad_P_dVi = 0;
821 PetscScalar d_grad_P_dVapp = 0;
822 PetscScalar d_grad_T_dTi = 0;
823
824 //--------------------------------
825 index[0] = zofs[z]+4*i+0;
826 index[1] = zofs[z]+4*i+1;
827 index[2] = zofs[z]+4*i+2;
828 index[3] = zofs[z]+4*i+3;
829 //--------------------------------
830 for(int j=0;j<pcell->nb_num;j++)
831 {
832 int nb = pcell->nb_array[j];
833 PetscScalar Vj = x[zofs[z]+4*nb+0]; //potential of nb node
834 PetscScalar nj = x[zofs[z]+4*nb+1]; //electron density of nb node
835 PetscScalar pj = x[zofs[z]+4*nb+2]; //hole density of nb node
836 PetscScalar Tj = x[zofs[z]+4*nb+3]; //lattice temperature of nb node
837 col[0] = zofs[z]+4*nb+0;
838 col[1] = zofs[z]+4*nb+1;
839 col[2] = zofs[z]+4*nb+2;
840 col[3] = zofs[z]+4*nb+3;
841 MatGetValues(*jtmp,4,index,4,col,A.m);
842 A=A/area;
843 if(j==0||j==pcell->nb_num-1)
844 {
845 PetscScalar Thick = pbc->Thick;
846 PetscScalar eps_ox = mt->eps0*pbc->eps;
847 PetscScalar eps = aux[i].eps;
848 PetscScalar s=eps_ox/eps/Thick;
849 d_grad_P_dVi += -0.5*pcell->ilen[j]*0.25*s*3/area;
850 d_grad_P_dVapp += 0.5*pcell->ilen[j]*s/area;
851 A.m[0]+= -0.5*pcell->ilen[j]*0.25*s/area;
852 PetscScalar h = pbc->heat_transfer;
853 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/area;
854 A.m[15] += -0.5*pcell->ilen[j]*0.25*h/area;
855 }
856 MatSetValues(*jac,4,index,4,col,A.m,INSERT_VALUES);
857
858 }
859 //fun(0) is the poisson's equation
860 //fun(1) is the continuous equation of electron
861 //fun(2) is the continuous equation of hole
862
863 MatGetValues(*jtmp,4,index,4,index,A.m);
864 A=A/area;
865 A.m[0] += d_grad_P_dVi;
866 A.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
867 A.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
868
869 A.m[15]+= d_grad_T_dTi; //dfun(3)/dT(i)
870 if(ODE_F.TimeDependent==true)
871 {
872 //second order
873 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
874 {
875 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
876 A.m[5] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
877 A.m[10] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
878 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
879 A.m[15] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
880 }
881 else //first order
882 {
883 A.m[5] += -1/ODE_F.dt;
884 A.m[10] += -1/ODE_F.dt;
885 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
886 A.m[15] += -aux[i].density*HeatCapacity/ODE_F.dt;
887 }
888 }
889 MatSetValues(*jac,4,index,4,index,A.m,INSERT_VALUES);
890 }
891
892
893 void SMCZone::F2E_mix_om_electrode_current(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
894 {
895 int bc_index = electrode[i];
896 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
897 PetscScalar e = mt->e;
898
899 //calculate the total current of ohmic electrode
900 PetscScalar current=0;
901 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
902 {
903 int node=bc[bc_index].psegment->node_array[j];
904 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
905 PetscScalar Vi = x[zofs[zone_index]+4*node+0]; //potential of node i
906 current += DeviceDepth*(f[zofs[zone_index]+4*node+1]-f[zofs[zone_index]+4*node+2]);
907 for(int k=0;k<pcell->nb_num;k++)
908 {
909 int nb = pcell->nb_array[k];
910 PetscScalar Vj = x[zofs[zone_index]+4*nb+0]; //potential of nb node
911 //displacement current
912 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
913 }
914 }
915 pbc->Set_Current_new(current);
916
917 }
918
919
920 void SMCZone::F2E_mix_stk_electrode_current(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
921 {
922 int bc_index = electrode[i];
923 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
924
925 PetscScalar current=0;
926 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
927 {
928 int node = bc[bc_index].psegment->node_array[j];
929 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
930 PetscScalar Vi = x[zofs[zone_index]+4*node+0]; //electron density of node i
931 PetscScalar ni = x[zofs[zone_index]+4*node+1]; //electron density of node i
932 PetscScalar pi = x[zofs[zone_index]+4*node+2]; //hole density of node i
933 PetscScalar Ti = x[zofs[zone_index]+4*node+3];
934 mt->mapping(&pzone->danode[node],&aux[node],ODE_F.clock);
935 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
936 PetscScalar Jsn = mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[node].affinity-deltaVB);
937 PetscScalar Jsp = mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[node].affinity+deltaVB);
938 current += -Jsn*0.5*pcell->ilen[0]*DeviceDepth;
939 current += -Jsp*0.5*pcell->ilen[0]*DeviceDepth;
940 current += -Jsn*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
941 current += -Jsp*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
942 for(int k=0;k<pcell->nb_num;k++)
943 {
944 int nb = pcell->nb_array[k];
945 PetscScalar Vj = x[zofs[zone_index]+4*nb+0]; //potential of nb node
946 //displacement current
947 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
948 }
949 }
950 pbc->Set_Current_new(current);
951
952 }
953
954
955 void SMCZone::F2E_mix_ins_electrode_current(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
956 {
957 int bc_index = electrode[i];
958 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
959 PetscScalar current=0;
960
961 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
962 {
963 int node = bc[bc_index].psegment->node_array[j];
964 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
965 PetscScalar Vi = x[zofs[zone_index]+4*node+0]; //electron density of node i
966 for(int k=0;k<pcell->nb_num;k++)
967 {
968 int nb = pcell->nb_array[k];
969 PetscScalar Vj = x[zofs[zone_index]+4*nb+0]; //potential of nb node
970 //displacement current
971 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
972 }
973 }
974 pbc->Set_Current_new(current);
975 }
976
977
978 void SMCZone::F2E_mix_om_electrode_Load(int bc_index, double & current, PetscScalar & pdI_pdV, Vec & pdI_pdw, Vec & pdF_pdV,PetscScalar *x, Mat *jtmp, ODE_Formula &ODE_F, vector<int> &zofs, DABC & bc, PetscScalar DeviceDepth)
979 {
980 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
981 pdI_pdV = 0;
982 PetscScalar e = mt->e;
983 PetscScalar kb = mt->kb;
984 VecZeroEntries(pdI_pdw);
985 VecZeroEntries(pdF_pdV);
986 PetscScalar * apdI_pdw;
987 PetscScalar * apdF_pdV;
988 VecGetArray(pdI_pdw,&apdI_pdw);
989 VecGetArray(pdF_pdV,&apdF_pdV);
990 current = pbc->Get_Current_new();
991 PetscScalar A1[4],A2[4];
992 PetscInt index[4],col[4];
993 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
994 {
995 int node=bc[bc_index].psegment->node_array[j];
996 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
997 PetscScalar Vi = x[zofs[zone_index]+4*node+0]; //potential of node i
998 PetscScalar dJdisp_dVi=0,dJdisp_dVr=0;
999 index[0] = zofs[zone_index]+4*node+0;
1000 index[1] = zofs[zone_index]+4*node+1;
1001 index[2] = zofs[zone_index]+4*node+2;
1002 index[3] = zofs[zone_index]+4*node+3;
1003 for(int k=0;k<pcell->nb_num;k++)
1004 {
1005 int nb = pcell->nb_array[k];
1006 PetscScalar Vj = x[zofs[zone_index]+4*nb+0]; //potential of nb node
1007 col[0] = zofs[zone_index]+4*nb+0;
1008 col[1] = zofs[zone_index]+4*nb+1;
1009 col[2] = zofs[zone_index]+4*nb+2;
1010 col[3] = zofs[zone_index]+4*nb+3;
1011
1012 MatGetValues(*jtmp,1,&index[1],4,col,A1);
1013 MatGetValues(*jtmp,1,&index[2],4,col,A2);
1014
1015 //for displacement current
1016 dJdisp_dVi += aux[node].eps/pcell->ilen[k]/ODE_F.dt*pcell->elen[k];
1017 dJdisp_dVr = -aux[node].eps/pcell->ilen[k]/ODE_F.dt*pcell->elen[k];
1018
1019 //
1020 apdI_pdw[col[0]] += (A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1021 apdI_pdw[col[1]] += (A1[1]-A2[1])*DeviceDepth;
1022 apdI_pdw[col[2]] += (A1[2]-A2[2])*DeviceDepth;
1023 apdI_pdw[col[3]] += (A1[3]-A2[3])*DeviceDepth;
1024 }
1025
1026 MatGetValues(*jtmp,1,&index[1],4,index,A1);
1027 MatGetValues(*jtmp,1,&index[2],4,index,A2);
1028
1029 apdI_pdw[index[0]] += (A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth;
1030 apdI_pdw[index[1]] += (A1[1]-A2[1])*DeviceDepth;
1031 apdI_pdw[index[2]] += (A1[2]-A2[2])*DeviceDepth;
1032 apdI_pdw[index[3]] += (A1[3]-A2[3])*DeviceDepth;
1033
1034 if(Fermi)
1035 {
1036 PetscScalar Vt = kb*fs[node].T/e;
1037 mt->mapping(&pzone->danode[node],&aux[node],ODE_F.clock);
1038 PetscScalar Nc = mt->band->Nc(fs[node].T);
1039 PetscScalar Nv = mt->band->Nv(fs[node].T);
1040 PetscScalar Vi = x[zofs[zone_index]+4*node+0]; //potential of node
1041 PetscScalar Ec = -(e*Vi + aux[node].affinity + mt->band->EgNarrowToEc(fs[node].T) + kb*fs[node].T*log(aux[node].Nc));
1042 PetscScalar Ev = -(e*Vi + aux[node].affinity - mt->band->EgNarrowToEv(fs[node].T) - kb*fs[node].T*log(aux[node].Nv)+ aux[node].Eg);
1043 PetscScalar phin = pbc->Vapp;
1044 PetscScalar phip = pbc->Vapp;
1045 PetscScalar etan = (-e*phin-Ec)/kb/fs[node].T;
1046 PetscScalar etap = (Ev+e*phip)/kb/fs[node].T;
1047 PetscScalar detan_dVi = 1.0/Vt;
1048 PetscScalar detap_dVi = -1.0/Vt;
1049 apdF_pdV[index[0]] = Nc*fermi_mhalf(etan)*detan_dVi - Nv*fermi_mhalf(etap)*detap_dVi;
1050 apdF_pdV[index[1]] = -Nc*fermi_mhalf(etan)*detan_dVi;
1051 apdF_pdV[index[2]] = -Nv*fermi_mhalf(etap)*detap_dVi;
1052 }
1053 else //Boltzmann
1054 apdF_pdV[index[0]] = 1.0;
1055 pdI_pdV = 0;
1056 }
1057 VecRestoreArray(pdI_pdw,&apdI_pdw);
1058 VecRestoreArray(pdF_pdV,&apdF_pdV);
1059 }
1060
1061
1062 void SMCZone::F2E_mix_stk_electrode_Load(int bc_index, double & current, PetscScalar & pdI_pdV, Vec & pdI_pdw, Vec & pdF_pdV,PetscScalar *x, Mat *jtmp, ODE_Formula &ODE_F, vector<int> &zofs, DABC & bc, PetscScalar DeviceDepth)
1063 {
1064 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
1065 pdI_pdV = 0;
1066 PetscScalar e = mt->e;
1067 VecZeroEntries(pdI_pdw);
1068 VecZeroEntries(pdF_pdV);
1069 PetscScalar * apdI_pdw;
1070 PetscScalar * apdF_pdV;
1071 VecGetArray(pdI_pdw,&apdI_pdw);
1072 VecGetArray(pdF_pdV,&apdF_pdV);
1073 current = pbc->Get_Current_new();
1074
1075 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1076 {
1077 int node = bc[bc_index].psegment->node_array[j];
1078 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1079 PetscScalar VB = pbc->WorkFunction-aux[node].affinity;
1080 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
1081 PetscScalar ni = x[zofs[zone_index]+4*node+1]; //electron density of node i
1082 PetscScalar pi = x[zofs[zone_index]+4*node+2]; //hole density of node i
1083 PetscScalar dJ_dni = -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[0]
1084 -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
1085 PetscScalar dJ_dpi = -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[0]
1086 -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
1087
1088 //for displacement current
1089 for(int k=0;k<pcell->nb_num;k++)
1090 {
1091 int nb = pcell->nb_array[k];
1092 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1093 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1094 apdI_pdw[zofs[zone_index]+4*node+0] += dI_dVi;
1095 apdI_pdw[zofs[zone_index]+4*nb+0] += dI_dVr;
1096 }
1097 apdI_pdw[zofs[zone_index]+4*node+1] = DeviceDepth*dJ_dni;
1098 apdI_pdw[zofs[zone_index]+4*node+2] = DeviceDepth*dJ_dpi;
1099
1100 apdF_pdV[zofs[zone_index]+4*node+0] = 1.0;
1101 pdI_pdV = 0;
1102 }
1103
1104 VecRestoreArray(pdI_pdw,&apdI_pdw);
1105 VecRestoreArray(pdF_pdV,&apdF_pdV);
1106 }
1107
1108
1109 void SMCZone::F2E_mix_ins_electrode_Load(int bc_index, double & current, PetscScalar & pdI_pdV, Vec & pdI_pdw, Vec & pdF_pdV,PetscScalar *x, Mat *jtmp, ODE_Formula &ODE_F, vector<int> &zofs, DABC & bc, PetscScalar DeviceDepth)
1110 {
1111 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
1112 pdI_pdV = 0;
1113 PetscScalar e = mt->e;
1114 VecZeroEntries(pdI_pdw);
1115 VecZeroEntries(pdF_pdV);
1116 PetscScalar * apdI_pdw;
1117 PetscScalar * apdF_pdV;
1118 VecGetArray(pdI_pdw,&apdI_pdw);
1119 VecGetArray(pdF_pdV,&apdF_pdV);
1120 current = pbc->Get_Current_new();
1121
1122 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1123 {
1124 int node = bc[bc_index].psegment->node_array[j];
1125 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1126
1127 //for displacement current
1128 for(int k=0;k<pcell->nb_num;k++)
1129 {
1130 int nb = pcell->nb_array[k];
1131 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1132 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1133 apdI_pdw[zofs[zone_index]+4*node+0] += dI_dVi;
1134 apdI_pdw[zofs[zone_index]+4*nb+0] += dI_dVr;
1135 if(k==0||k==pcell->nb_num-1)
1136 {
1137 PetscScalar Thick = pbc->Thick;
1138 PetscScalar eps_ox = mt->eps0*pbc->eps;
1139 PetscScalar s=eps_ox/Thick;
1140 apdF_pdV[zofs[zone_index]+4*node+0] += 0.5*pcell->ilen[k]*s/pcell->area;
1141 }
1142 }
1143 pdI_pdV = 0;
1144 }
1145
1146 VecRestoreArray(pdI_pdw,&apdI_pdw);
1147 VecRestoreArray(pdF_pdV,&apdF_pdV);
1148 }
1149
General-purpose Semiconductor Simulator