?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix: fixes Ec/Ev in all solver.
[gss-tcad.git] / src / solver / ebm3e / semiequ3e.cc
blob9c984a5e23a5af4346b2cdc3d08362e90c310ac9
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 15, 2007 */
14 /* */
15 /* Gong Ding */
16 /* gdiso@ustc.edu */
17 /* NINT, No.69 P.O.Box, Xi'an City, China */
18 /* */
19 /*****************************************************************************/
20
21 #include "mathfunc.h"
22 #include "flux3e.h"
23 #include "zonedata.h"
24 #include "vec6.h"
25
26 #define _HMOB_
27 #define _APPROXIMATE_
28
29
30 void SMCZone::F3E_Tri_ddm(Tri *ptri,PetscScalar *x,PetscScalar *f, vector<int> &zofs)
31 {
32 PetscScalar grad_P,sn,sp,kapa,grad_T,H=0,S;
33 PetscScalar Epn,Epp,Etn,Etp;
34 int A = ptri->node[0];
35 int B = ptri->node[1];
36 int C = ptri->node[2];
37 PetscScalar xa = pzone->danode[ptri->node[0]].x;
38 PetscScalar xb = pzone->danode[ptri->node[1]].x;
39 PetscScalar xc = pzone->danode[ptri->node[2]].x;
40 PetscScalar ya = pzone->danode[ptri->node[0]].y;
41 PetscScalar yb = pzone->danode[ptri->node[1]].y;
42 PetscScalar yc = pzone->danode[ptri->node[2]].y;
43 PetscScalar tri_area = ptri->area;
44 PetscScalar Sax = (xc-xb)/ptri->edge_len[0];
45 PetscScalar Say = (yc-yb)/ptri->edge_len[0];
46 PetscScalar Sbx = (xa-xc)/ptri->edge_len[1];
47 PetscScalar Sby = (ya-yc)/ptri->edge_len[1];
48 PetscScalar Scx = (xb-xa)/ptri->edge_len[2];
49 PetscScalar Scy = (yb-ya)/ptri->edge_len[2];
50 PetscScalar Wab = ptri->d[1]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
51 PetscScalar Wac = ptri->d[2]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[1])*sin(ptri->angle[1]));
52 PetscScalar Wbc = ptri->d[2]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[0])*sin(ptri->angle[0]));
53 PetscScalar Wba = ptri->d[0]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
54 PetscScalar Wca = ptri->d[0]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[1])*sin(ptri->angle[1]));
55 PetscScalar Wcb = ptri->d[1]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[0])*sin(ptri->angle[0]));
56 PetscScalar Ma = -cos(ptri->angle[0]);
57 PetscScalar Mb = -cos(ptri->angle[1]);
58 PetscScalar Mc = -cos(ptri->angle[2]);
59
60 PetscScalar e = mt->e;
61 PetscScalar kb = mt->kb;
62 PetscScalar Tl = (fs[A].T+fs[B].T+fs[C].T)/3.0;
63
64 PetscScalar Va = x[zofs[zone_index]+6*A+0]; //potential of node A
65 PetscScalar na = x[zofs[zone_index]+6*A+1]; //electron density of node A
66 PetscScalar pa = x[zofs[zone_index]+6*A+2]; //hole density of node A
67 PetscScalar Ta = x[zofs[zone_index]+6*A+3]; //Lattice Temperature of node A
68 PetscScalar Tna = x[zofs[zone_index]+6*A+4]/na; //elec temperature of node A
69 PetscScalar Tpa = x[zofs[zone_index]+6*A+5]/pa; //hole temperature of node A
70 mt->mapping(&pzone->danode[A],&aux[A],0);
71 PetscScalar nia = mt->band->nie(Ta);
72 PetscScalar Ra = mt->band->Recomb(pa,na,Ta);
73 PetscScalar Ra_SHR = mt->band->R_SHR(pa,na,Ta);
74 PetscScalar Ra_AUG = mt->band->R_Auger(pa,na,Ta);
75 PetscScalar Ra_DIR = mt->band->R_Direct(pa,na,Ta);
76 #ifdef _APPROXIMATE_
77 PetscScalar Nca = mt->band->Nc(fs[A].T);
78 PetscScalar Nva = mt->band->Nv(fs[A].T);
79 PetscScalar Ega = mt->band->Eg(fs[A].T);
80 PetscScalar Eca = -(e*Va + aux[A].affinity + mt->band->EgNarrowToEc(fs[A].T) + kb*fs[A].T*(log(Nca)-1.5*log(fs[A].T)));
81 PetscScalar Eva = -(e*Va + aux[A].affinity - mt->band->EgNarrowToEv(fs[A].T) - kb*fs[A].T*(log(Nva)-1.5*log(fs[A].T)) + Ega);
82 #else
83 PetscScalar Nca = mt->band->Nc(Ta);
84 PetscScalar Nva = mt->band->Nv(Ta);
85 PetscScalar Ega = mt->band->Eg(Ta);
86 PetscScalar Eca = -(e*Va + aux[A].affinity + mt->band->EgNarrowToEc(Ta) + kb*Ta*(log(Nca)-1.5*log(Ta)));
87 PetscScalar Eva = -(e*Va + aux[A].affinity - mt->band->EgNarrowToEv(Ta) - kb*Ta*(log(Nva)-1.5*log(Ta)) + Ega);
88 #endif
89 PetscScalar tao_ena = mt->band->ElecEnergyRelaxTime(Tna,Ta);
90 PetscScalar tao_epa = mt->band->HoleEnergyRelaxTime(Tpa,Ta);
91
92 PetscScalar Vb = x[zofs[zone_index]+6*B+0]; //potential of node B
93 PetscScalar nb = x[zofs[zone_index]+6*B+1]; //electron density of node B
94 PetscScalar pb = x[zofs[zone_index]+6*B+2]; //hole density of node B
95 PetscScalar Tb = x[zofs[zone_index]+6*B+3]; //Temperature of node B
96 PetscScalar Tnb = x[zofs[zone_index]+6*B+4]/nb; //elec temperature of node B
97 PetscScalar Tpb = x[zofs[zone_index]+6*B+5]/pb; //hole temperature of node B
98 mt->mapping(&pzone->danode[B],&aux[B],0);
99 PetscScalar nib = mt->band->nie(Tb);
100 PetscScalar Rb = mt->band->Recomb(pb,nb,Tb);
101 PetscScalar Rb_SHR = mt->band->R_SHR(pb,nb,Tb);
102 PetscScalar Rb_AUG = mt->band->R_Auger(pb,nb,Tb);
103 PetscScalar Rb_DIR = mt->band->R_Direct(pb,nb,Tb);
104 #ifdef _APPROXIMATE_
105 PetscScalar Ncb = mt->band->Nc(fs[B].T);
106 PetscScalar Nvb = mt->band->Nv(fs[B].T);
107 PetscScalar Egb = mt->band->Eg(fs[B].T);
108 PetscScalar Ecb = -(e*Vb + aux[B].affinity + mt->band->EgNarrowToEc(fs[B].T) + kb*fs[B].T*(log(Ncb)-1.5*log(fs[B].T)));
109 PetscScalar Evb = -(e*Vb + aux[B].affinity - mt->band->EgNarrowToEv(fs[B].T) - kb*fs[B].T*(log(Nvb)-1.5*log(fs[B].T))+ Egb);
110 #else
111 PetscScalar Ncb = mt->band->Nc(Tb);
112 PetscScalar Nvb = mt->band->Nv(Tb);
113 PetscScalar Egb = mt->band->Eg(Tb);
114 PetscScalar Ecb = -(e*Vb + aux[B].affinity + mt->band->EgNarrowToEc(Tb) + kb*Tb*(log(Ncb)-1.5*log(Tb)));
115 PetscScalar Evb = -(e*Vb + aux[B].affinity - mt->band->EgNarrowToEv(Tb) - kb*Tb*(log(Nvb)-1.5*log(Tb))+ Egb);
116 #endif
117 PetscScalar tao_enb = mt->band->ElecEnergyRelaxTime(Tnb,Tb);
118 PetscScalar tao_epb = mt->band->HoleEnergyRelaxTime(Tpb,Tb);
119
120 PetscScalar Vc = x[zofs[zone_index]+6*C+0]; //potential of node C
121 PetscScalar nc = x[zofs[zone_index]+6*C+1]; //electron density of node C
122 PetscScalar pc = x[zofs[zone_index]+6*C+2]; //hole density of node C
123 PetscScalar Tc = x[zofs[zone_index]+6*C+3]; //Temperature of node C
124 PetscScalar Tnc = x[zofs[zone_index]+6*C+4]/nc; //Temperature of node C
125 PetscScalar Tpc = x[zofs[zone_index]+6*C+5]/pc; //Temperature of node C
126 mt->mapping(&pzone->danode[C],&aux[C],0);
127 PetscScalar nic = mt->band->nie(Tc);
128 PetscScalar Rc = mt->band->Recomb(pc,nc,Tc);
129 PetscScalar Rc_SHR = mt->band->R_SHR(pc,nc,Tc);
130 PetscScalar Rc_AUG = mt->band->R_Auger(pc,nc,Tc);
131 PetscScalar Rc_DIR = mt->band->R_Direct(pc,nc,Tc);
132 #ifdef _APPROXIMATE_
133 PetscScalar Ncc = mt->band->Nc(fs[C].T);
134 PetscScalar Nvc = mt->band->Nv(fs[C].T);
135 PetscScalar Egc = mt->band->Eg(fs[C].T);
136 PetscScalar Ecc = -(e*Vc + aux[C].affinity + mt->band->EgNarrowToEc(fs[C].T) + kb*fs[C].T*(log(Ncc)-1.5*log(fs[C].T)));
137 PetscScalar Evc = -(e*Vc + aux[C].affinity - mt->band->EgNarrowToEv(fs[C].T) - kb*fs[C].T*(log(Nvc)-1.5*log(fs[C].T))+ Egc);
138 #else
139 PetscScalar Ncc = mt->band->Nc(Tc);
140 PetscScalar Nvc = mt->band->Nv(Tc);
141 PetscScalar Egc = mt->band->Eg(Tc);
142 PetscScalar Ecc = -(e*Vc + aux[C].affinity + mt->band->EgNarrowToEc(Tc) + kb*Tc*(log(Ncc)-1.5*log(Tc)));
143 PetscScalar Evc = -(e*Vc + aux[C].affinity - mt->band->EgNarrowToEv(Tc) - kb*Tc*(log(Nvc)-1.5*log(Tc))+ Egc);
144 #endif
145 PetscScalar tao_enc = mt->band->ElecEnergyRelaxTime(Tnc,Tc);
146 PetscScalar tao_epc = mt->band->HoleEnergyRelaxTime(Tpc,Tc);
147
148 PetscScalar Ex = -((yb-yc)*Va + (yc-ya)*Vb +(ya-yb)*Vc)/(2*tri_area);
149 PetscScalar Ey = -((xc-xb)*Va + (xa-xc)*Vb +(xb-xa)*Vc)/(2*tri_area);
150 PetscScalar E = sqrt(Ex*Ex+Ey*Ey+1e-20);
151 PetscScalar Eg = mt->band->Eg((Ta+Tb+Tc)/3.0);
152
153 PetscScalar Ga=0,Gb=0,Gc=0,G=0;
154 PetscScalar IIn=0,IIp=0;
155 PetscScalar mun,mup;
156 PetscScalar Epnc=0,Etnc=0;
157 PetscScalar Eppc=0,Etpc=0;
158 PetscScalar Epna=0,Etna=0;
159 PetscScalar Eppa=0,Etpa=0;
160 PetscScalar Epnb=0,Etnb=0;
161 PetscScalar Eppb=0,Etpb=0;
162 PetscScalar Jna_norm,Jpa_norm;
163 PetscScalar Jnb_norm,Jpb_norm;
164 PetscScalar Jnc_norm,Jpc_norm;
165
166 PetscScalar Jnc = In(kb, e, -Eca/e, -Ecb/e, na, nb, Tna,Tnb, ptri->edge_len[2]);
167 PetscScalar Jpc = Ip(kb, e, -Eva/e, -Evb/e, pa, pb, Tpa,Tpb, ptri->edge_len[2]);
168 PetscScalar Jna = In(kb, e, -Ecb/e, -Ecc/e, nb, nc, Tnb,Tnc, ptri->edge_len[0]);
169 PetscScalar Jpa = Ip(kb, e, -Evb/e, -Evc/e, pb, pc, Tpb,Tpc, ptri->edge_len[0]);
170 PetscScalar Jnb = In(kb, e, -Ecc/e, -Eca/e, nc, na, Tnc,Tna, ptri->edge_len[1]);
171 PetscScalar Jpb = Ip(kb, e, -Evc/e, -Eva/e, pc, pa, Tpc,Tpa, ptri->edge_len[1]);
172 PetscScalar Jn_scale = dmax(dmax(fabs(Jna),fabs(Jnb)),dmax(fabs(Jnc),nia*nia));
173 PetscScalar Jp_scale = dmax(dmax(fabs(Jpa),fabs(Jpb)),dmax(fabs(Jpc),nia*nia));
174 Jnc /= Jn_scale;
175 Jpc /= Jp_scale;
176 Jna /= Jn_scale;
177 Jpa /= Jp_scale;
178 Jnb /= Jn_scale;
179 Jpb /= Jp_scale;
180
181 //flux along A-B
182 if(HighFieldMobility)
183 {
184 if(EJModel || IIType==EdotJ || IIType==TempII)
185 {
186 PetscScalar Jncx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jnc
187 + (Wca*Sax - Wca*Mb*Scx)*Jna
188 + (Wcb*Sbx - Wcb*Ma*Scx)*Jnb;
189 PetscScalar Jncy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jnc
190 + (Wca*Say - Wca*Mb*Scy)*Jna
191 + (Wcb*Sby - Wcb*Ma*Scy)*Jnb;
192 PetscScalar Jpcx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jpc
193 + (Wca*Sax - Wca*Mb*Scx)*Jpa
194 + (Wcb*Sbx - Wcb*Ma*Scx)*Jpb;
195 PetscScalar Jpcy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jpc
196 + (Wca*Say - Wca*Mb*Scy)*Jpa
197 + (Wcb*Sby - Wcb*Ma*Scy)*Jpb;
198 Jnc_norm = sqrt(Jncx*Jncx + Jncy*Jncy + 1e-100);
199 Jpc_norm = sqrt(Jpcx*Jpcx + Jpcy*Jpcy + 1e-100);
200 Epnc = (Ex*Jncx + Ey*Jncy)/Jnc_norm;
201 Etnc = (Ex*Jncy - Ey*Jncx)/Jnc_norm;
202 Eppc = (Ex*Jpcx + Ey*Jpcy)/Jpc_norm;
203 Etpc = (Ex*Jpcy - Ey*Jpcx)/Jpc_norm;
204 }
205 else
206 {
207 Epnc = fabs(Eca-Ecb)/e/ptri->edge_len[2]; //parallel electrical field for electron
208 Eppc = fabs(Eva-Evb)/e/ptri->edge_len[2]; //parallel electrical field for hole
209 //transvers electrical field for electron and hole
210 Etnc = fabs(Eca + ptri->edge_len[1]*cos(ptri->angle[0])*(Ecb-Eca)/ptri->edge_len[2] - Ecc)
211 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
212 Etpc = fabs(Eva + ptri->edge_len[1]*cos(ptri->angle[0])*(Evb-Eva)/ptri->edge_len[2] - Evc)
213 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
214 }
215
216 mt->mapping(&pzone->danode[A],&aux[A],0);
217 aux[A].mun = mt->mob->ElecMob(pa,na,Ta,dmax(0,Epnc),fabs(Etnc),1);
218 aux[A].mup = mt->mob->HoleMob(pa,na,Ta,dmax(0,Eppc),fabs(Etpc),1);
219
220 mt->mapping(&pzone->danode[B],&aux[B],0);
221 aux[B].mun = mt->mob->ElecMob(pb,nb,Tb,dmax(0,Epnc),fabs(Etnc),1);
222 aux[B].mup = mt->mob->HoleMob(pb,nb,Tb,dmax(0,Eppc),fabs(Etpc),1);
223 }
224
225 mun = 0.5*(aux[A].mun+aux[B].mun);
226 mup = 0.5*(aux[A].mup+aux[B].mup);
227 if(ImpactIonization)
228 {
229 if(IIType==EdotJ)
230 {
231 IIn = mt->gen->ElecGenRate(0.5*(Ta+Tb),dmax(0,Epnc),Eg);
232 IIp = mt->gen->HoleGenRate(0.5*(Ta+Tb),dmax(0,Eppc),Eg);
233 Gc = IIn*mun*Jnc_norm*Jn_scale + IIp*mup*Jpc_norm*Jp_scale;
234 }
235 if(IIType==TempII)
236 {
237 IIn = mt->gen->ElecGenRateEBM(0.5*(Tna+Tnb),0.5*(Ta+Tb),Eg);
238 IIp = mt->gen->HoleGenRateEBM(0.5*(Tpa+Tpb),0.5*(Ta+Tb),Eg);
239 Gc = IIn*mun*Jnc_norm*Jn_scale + IIp*mup*Jpc_norm*Jp_scale;
240 }
241 if(IIType==ESide )
242 {
243 IIn = mt->gen->ElecGenRate(0.5*(Ta+Tb),fabs(Va-Vb)/ptri->edge_len[2],Eg);
244 IIp = mt->gen->HoleGenRate(0.5*(Ta+Tb),fabs(Va-Vb)/ptri->edge_len[2],Eg);
245 Gc = IIn*mun*fabs(Jnc)*Jn_scale + IIp*mup*fabs(Jpc)*Jp_scale;
246 }
247 }
248 grad_P = (Vb-Va)/ptri->edge_len[2];
249 kapa = mt->thermal->HeatConduction(0.5*(Ta+Tb));
250 grad_T = kapa*(Tb-Ta)/ptri->edge_len[2];
251
252 sn = Sn(kb, e, -Eca/e, -Ecb/e, na, nb, Tna,Tnb, ptri->edge_len[2]);
253 sp = Sp(kb, e, -Eva/e, -Evb/e, pa, pb, Tpa,Tpb, ptri->edge_len[2]);
254
255 f[zofs[zone_index]+6*A+0] += grad_P*ptri->d[2];
256 f[zofs[zone_index]+6*A+1] += mun*Jnc*Jn_scale*ptri->d[2] + Gc*0.25*ptri->d[2]*ptri->edge_len[2];
257 f[zofs[zone_index]+6*A+2] += - mup*Jpc*Jp_scale*ptri->d[2] + Gc*0.25*ptri->d[2]*ptri->edge_len[2];
258 f[zofs[zone_index]+6*A+3] += grad_T*ptri->d[2];
259 f[zofs[zone_index]+6*A+4] += - mun*sn*ptri->d[2] + 0.5*(Va-Vb)*mun*Jnc*Jn_scale*ptri->d[2];
260 f[zofs[zone_index]+6*A+5] += - mup*sp*ptri->d[2] + 0.5*(Va-Vb)*mup*Jpc*Jp_scale*ptri->d[2];
261
262 f[zofs[zone_index]+6*B+0] += - grad_P*ptri->d[2];
263 f[zofs[zone_index]+6*B+1] += - mun*Jnc*Jn_scale*ptri->d[2] + Gc*0.25*ptri->d[2]*ptri->edge_len[2];
264 f[zofs[zone_index]+6*B+2] += mup*Jpc*Jp_scale*ptri->d[2] + Gc*0.25*ptri->d[2]*ptri->edge_len[2];
265 f[zofs[zone_index]+6*B+3] += - grad_T*ptri->d[2];
266 f[zofs[zone_index]+6*B+4] += mun*sn*ptri->d[2] + 0.5*(Va-Vb)*mun*Jnc*Jn_scale*ptri->d[2];
267 f[zofs[zone_index]+6*B+5] += mup*sp*ptri->d[2] + 0.5*(Va-Vb)*mup*Jpc*Jp_scale*ptri->d[2];
268
269 //flux along B-C
270 if(HighFieldMobility)
271 {
272 if(EJModel || IIType==EdotJ || IIType==TempII)
273 {
274 PetscScalar Jnax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jna
275 + (Wab*Sbx - Wab*Mc*Sax)*Jnb
276 + (Wac*Scx - Wac*Mb*Sax)*Jnc;
277 PetscScalar Jnay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jna
278 + (Wab*Sby - Wab*Mc*Say)*Jnb
279 + (Wac*Scy - Wac*Mb*Say)*Jnc;
280 PetscScalar Jpax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jpa
281 + (Wab*Sbx - Wab*Mc*Sax)*Jpb
282 + (Wac*Scx - Wac*Mb*Sax)*Jpc;
283 PetscScalar Jpay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jpa
284 + (Wab*Sby - Wab*Mc*Say)*Jpb
285 + (Wac*Scy - Wac*Mb*Say)*Jpc;
286 Jna_norm = sqrt(Jnax*Jnax + Jnay*Jnay + 1e-100);
287 Jpa_norm = sqrt(Jpax*Jpax + Jpay*Jpay + 1e-100);
288 Epna = (Ex*Jnax + Ey*Jnay)/Jna_norm;
289 Etna = (Ex*Jnay - Ey*Jnax)/Jna_norm;
290 Eppa = (Ex*Jpax + Ey*Jpay)/Jpa_norm;
291 Etpa = (Ex*Jpay - Ey*Jpax)/Jpa_norm;
292 }
293 else
294 {
295 Epna = fabs(Ecb-Ecc)/e/ptri->edge_len[0]; //parallel electrical field for electron
296 Eppa = fabs(Evb-Evc)/e/ptri->edge_len[0]; //parallel electrical field for hole
297 //transvers electrical field for electron and hole
298 Etna = fabs(Ecb + ptri->edge_len[2]*cos(ptri->angle[1])*(Ecc-Ecb)/ptri->edge_len[0] - Eca)
299 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
300 Etpa = fabs(Evb + ptri->edge_len[2]*cos(ptri->angle[1])*(Evc-Evb)/ptri->edge_len[0] - Eva)
301 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
302 }
303 mt->mapping(&pzone->danode[B],&aux[B],0);
304 aux[B].mun = mt->mob->ElecMob(pb,nb,Tb,dmax(0,Epna),fabs(Etna),1);
305 aux[B].mup = mt->mob->HoleMob(pb,nb,Tb,dmax(0,Eppa),fabs(Etpa),1);
306
307 mt->mapping(&pzone->danode[C],&aux[C],0);
308 aux[C].mun = mt->mob->ElecMob(pc,nc,Tc,dmax(0,Epna),fabs(Etna),1);
309 aux[C].mup = mt->mob->HoleMob(pc,nc,Tc,dmax(0,Eppa),fabs(Etpa),1);
310 }
311
312 mun = 0.5*(aux[B].mun+aux[C].mun);
313 mup = 0.5*(aux[B].mup+aux[C].mup);
314
315 if(ImpactIonization)
316 {
317 if(IIType==EdotJ)
318 {
319 IIn = mt->gen->ElecGenRate(0.5*(Tb+Tc),dmax(0,Epna),Eg);
320 IIp = mt->gen->HoleGenRate(0.5*(Tb+Tc),dmax(0,Eppa),Eg);
321 Ga = IIn*mun*Jna_norm*Jn_scale + IIp*mup*Jpa_norm*Jp_scale;
322 }
323 if(IIType==TempII)
324 {
325 IIn = mt->gen->ElecGenRateEBM(0.5*(Tnb+Tnc),0.5*(Tb+Tc),Eg);
326 IIp = mt->gen->HoleGenRateEBM(0.5*(Tpb+Tpc),0.5*(Tb+Tc),Eg);
327 Ga = IIn*mun*Jna_norm*Jn_scale + IIp*mup*Jpa_norm*Jp_scale;
328 }
329 if(IIType==ESide)
330 {
331 IIn = mt->gen->ElecGenRate(0.5*(Tb+Tc),fabs(Vb-Vc)/ptri->edge_len[0],Eg);
332 IIp = mt->gen->HoleGenRate(0.5*(Tb+Tc),fabs(Vb-Vc)/ptri->edge_len[0],Eg);
333 Ga = IIn*mun*fabs(Jna)*Jn_scale + IIp*mup*fabs(Jpa)*Jp_scale;
334 }
335 }
336
337 grad_P = (Vc-Vb)/ptri->edge_len[0];
338 kapa = mt->thermal->HeatConduction(0.5*(Tb+Tc));
339 grad_T = kapa*(Tc-Tb)/ptri->edge_len[0];
340
341 sn = Sn(kb, e, -Ecb/e, -Ecc/e, nb, nc, Tnb,Tnc, ptri->edge_len[0]);
342 sp = Sp(kb, e, -Evb/e, -Evc/e, pb, pc, Tpb,Tpc, ptri->edge_len[0]);
343
344 f[zofs[zone_index]+6*B+0] += grad_P*ptri->d[0];
345 f[zofs[zone_index]+6*B+1] += mun*Jna*Jn_scale*ptri->d[0] + Ga*0.25*ptri->d[0]*ptri->edge_len[0];
346 f[zofs[zone_index]+6*B+2] += - mup*Jpa*Jp_scale*ptri->d[0] + Ga*0.25*ptri->d[0]*ptri->edge_len[0];
347 f[zofs[zone_index]+6*B+3] += grad_T*ptri->d[0];
348 f[zofs[zone_index]+6*B+4] += - mun*sn*ptri->d[0] + 0.5*(Vb-Vc)*mun*Jna*Jn_scale*ptri->d[0];
349 f[zofs[zone_index]+6*B+5] += - mup*sp*ptri->d[0] + 0.5*(Vb-Vc)*mup*Jpa*Jp_scale*ptri->d[0];
350
351 f[zofs[zone_index]+6*C+0] += - grad_P*ptri->d[0];
352 f[zofs[zone_index]+6*C+1] += - mun*Jna*Jn_scale*ptri->d[0] + Ga*0.25*ptri->d[0]*ptri->edge_len[0];
353 f[zofs[zone_index]+6*C+2] += mup*Jpa*Jp_scale*ptri->d[0] + Ga*0.25*ptri->d[0]*ptri->edge_len[0];
354 f[zofs[zone_index]+6*C+3] += - grad_T*ptri->d[0];
355 f[zofs[zone_index]+6*C+4] += mun*sn*ptri->d[0] + 0.5*(Vb-Vc)*mun*Jna*Jn_scale*ptri->d[0];
356 f[zofs[zone_index]+6*C+5] += mup*sp*ptri->d[0] + 0.5*(Vb-Vc)*mup*Jpa*Jp_scale*ptri->d[0];
357
358
359 //flux along C-A
360 if(HighFieldMobility)
361 {
362 if(EJModel || IIType==EdotJ || IIType==TempII)
363 {
364 PetscScalar Jnbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jnb
365 + (Wbc*Scx - Wbc*Ma*Sbx)*Jnc
366 + (Wba*Sax - Wba*Mc*Sbx)*Jna;
367 PetscScalar Jnby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jnb
368 + (Wbc*Scy - Wbc*Ma*Sby)*Jnc
369 + (Wba*Say - Wba*Mc*Sby)*Jna;
370 PetscScalar Jpbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jpb
371 + (Wbc*Scx - Wbc*Ma*Sbx)*Jpc
372 + (Wba*Sax - Wba*Mc*Sbx)*Jpa;
373 PetscScalar Jpby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jpb
374 + (Wbc*Scy - Wbc*Ma*Sby)*Jpc
375 + (Wba*Say - Wba*Mc*Sby)*Jpa;
376 Jnb_norm = sqrt(Jnbx*Jnbx + Jnby*Jnby + 1e-100);
377 Jpb_norm = sqrt(Jpbx*Jpbx + Jpby*Jpby + 1e-100);
378 Epnb = (Ex*Jnbx + Ey*Jnby)/Jnb_norm;
379 Etnb = (Ex*Jnby - Ey*Jnbx)/Jnb_norm;
380 Eppb = (Ex*Jpbx + Ey*Jpby)/Jpb_norm;
381 Etpb = (Ex*Jpby - Ey*Jpbx)/Jpb_norm;
382
383 }
384 else
385 {
386 Epnb = fabs(Ecc-Eca)/e/ptri->edge_len[1]; //parallel electrical field for electron
387 Eppb = fabs(Evc-Eva)/e/ptri->edge_len[1]; //parallel electrical field for hole
388 //transvers electrical field for electron and hole
389 Etnb = fabs(Ecc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eca-Ecc)/ptri->edge_len[1] - Ecb)
390 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
391 Etpb = fabs(Evc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eva-Evc)/ptri->edge_len[1] - Evb)
392 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
393 }
394 mt->mapping(&pzone->danode[C],&aux[C],0);
395 aux[C].mun = mt->mob->ElecMob(pc,nc,Tc,dmax(0,Epnb),fabs(Etnb),1);
396 aux[C].mup = mt->mob->HoleMob(pc,nc,Tc,dmax(0,Eppb),fabs(Etpb),1);
397
398 mt->mapping(&pzone->danode[A],&aux[A],0);
399 aux[A].mun = mt->mob->ElecMob(pa,na,Ta,dmax(0,Epnb),fabs(Etnb),1);
400 aux[A].mup = mt->mob->HoleMob(pa,na,Ta,dmax(0,Eppb),fabs(Etpb),1);
401 }
402
403 mun = 0.5*(aux[C].mun+aux[A].mun);
404 mup = 0.5*(aux[C].mup+aux[A].mup);
405
406 if(ImpactIonization)
407 {
408 if(IIType==EdotJ)
409 {
410 IIn = mt->gen->ElecGenRate(0.5*(Tc+Ta),dmax(0,Epnb),Eg);
411 IIp = mt->gen->HoleGenRate(0.5*(Tc+Ta),dmax(0,Eppb),Eg);
412 Gb = IIn*mun*Jnb_norm*Jn_scale + IIp*mup*Jpb_norm*Jp_scale;
413 }
414 if(IIType==TempII)
415 {
416 IIn = mt->gen->ElecGenRateEBM(0.5*(Tnc+Tna),0.5*(Tc+Ta),Eg);
417 IIp = mt->gen->HoleGenRateEBM(0.5*(Tpc+Tpa),0.5*(Tc+Ta),Eg);
418 Gb = IIn*mun*Jnb_norm*Jn_scale + IIp*mup*Jpb_norm*Jp_scale;
419 }
420 if(IIType==ESide)
421 {
422 IIn = mt->gen->ElecGenRate(0.5*(Tc+Ta),fabs(Vc-Va)/ptri->edge_len[1],Eg);
423 IIp = mt->gen->HoleGenRate(0.5*(Tc+Ta),fabs(Vc-Va)/ptri->edge_len[1],Eg);
424 Gb = IIn*mun*fabs(Jnb)*Jn_scale + IIp*mup*fabs(Jpb)*Jp_scale;
425 }
426 }
427
428 grad_P = (Va-Vc)/ptri->edge_len[1];
429 kapa = mt->thermal->HeatConduction(0.5*(Tc+Ta));
430 grad_T = kapa*(Ta-Tc)/ptri->edge_len[1];
431
432 sn = Sn(kb, e, -Ecc/e, -Eca/e, nc, na, Tnc,Tna, ptri->edge_len[1]);
433 sp = Sp(kb, e, -Evc/e, -Eva/e, pc, pa, Tpc,Tpa, ptri->edge_len[1]);
434
435 f[zofs[zone_index]+6*C+0] += grad_P*ptri->d[1];
436 f[zofs[zone_index]+6*C+1] += mun*Jnb*Jn_scale*ptri->d[1] + Gb*0.25*ptri->d[1]*ptri->edge_len[1];
437 f[zofs[zone_index]+6*C+2] += - mup*Jpb*Jp_scale*ptri->d[1] + Gb*0.25*ptri->d[1]*ptri->edge_len[1];
438 f[zofs[zone_index]+6*C+3] += grad_T*ptri->d[1];
439 f[zofs[zone_index]+6*C+4] += - mun*sn*ptri->d[1] + 0.5*(Vc-Va)*mun*Jnb*Jn_scale*ptri->d[1];
440 f[zofs[zone_index]+6*C+5] += - mup*sp*ptri->d[1] + 0.5*(Vc-Va)*mup*Jpb*Jp_scale*ptri->d[1];
441
442 f[zofs[zone_index]+6*A+0] += - grad_P*ptri->d[1];
443 f[zofs[zone_index]+6*A+1] += - mun*Jnb*Jn_scale*ptri->d[1] + Gb*0.25*ptri->d[1]*ptri->edge_len[1];
444 f[zofs[zone_index]+6*A+2] += mup*Jpb*Jp_scale*ptri->d[1] + Gb*0.25*ptri->d[1]*ptri->edge_len[1];
445 f[zofs[zone_index]+6*A+3] += - grad_T*ptri->d[1];
446 f[zofs[zone_index]+6*A+4] += mun*sn*ptri->d[1] + 0.5*(Vc-Va)*mun*Jnb*Jn_scale*ptri->d[1];
447 f[zofs[zone_index]+6*A+5] += mup*sp*ptri->d[1] + 0.5*(Vc-Va)*mup*Jpb*Jp_scale*ptri->d[1];
448
449 //---------------------------------------------------------------------------
450 // G-R item and heat source for each node
451 //---------------------------------------------------------------------------
452 G=(Ga+Gb+Gc)/3.0;
453 S = 0.25*ptri->d[2]*ptri->edge_len[2] + 0.25*ptri->d[1]*ptri->edge_len[1];
454 PetscScalar Hna = (Ra_AUG-G)*(Ega+1.5*kb*Tpa) - 1.5*kb*Tna*(Ra_SHR+Ra_DIR-G) - 1.5*kb*na*(Tna-Ta)/tao_ena;
455 PetscScalar Hpa = (Ra_AUG-G)*(Ega+1.5*kb*Tna) - 1.5*kb*Tpa*(Ra_SHR+Ra_DIR-G) - 1.5*kb*pa*(Tpa-Ta)/tao_epa;
456 PetscScalar Ha = Ra_SHR*(Ega+1.5*kb*Tna+1.5*kb*Tpa) + 1.5*kb*na*(Tna-Ta)/tao_ena + 1.5*kb*pa*(Tpa-Ta)/tao_epa;
457 f[zofs[zone_index]+6*A+1] += (G-Ra)*S;
458 f[zofs[zone_index]+6*A+2] += (G-Ra)*S;
459 f[zofs[zone_index]+6*A+3] += Ha*S;
460 f[zofs[zone_index]+6*A+4] += Hna*S;
461 f[zofs[zone_index]+6*A+5] += Hpa*S;
462
463 S = 0.25*ptri->d[0]*ptri->edge_len[0] + 0.25*ptri->d[2]*ptri->edge_len[2];
464 PetscScalar Hnb = (Rb_AUG-G)*(Egb+1.5*kb*Tpb) - 1.5*kb*Tnb*(Rb_SHR+Rb_DIR-G) - 1.5*kb*nb*(Tnb-Tb)/tao_enb;
465 PetscScalar Hpb = (Rb_AUG-G)*(Egb+1.5*kb*Tnb) - 1.5*kb*Tpb*(Rb_SHR+Rb_DIR-G) - 1.5*kb*pb*(Tpb-Tb)/tao_epb;
466 PetscScalar Hb = Rb_SHR*(Egb+1.5*kb*Tnb+1.5*kb*Tpb) + 1.5*kb*nb*(Tnb-Tb)/tao_enb + 1.5*kb*pb*(Tpb-Tb)/tao_epb;
467 f[zofs[zone_index]+6*B+1] += (G-Rb)*S;
468 f[zofs[zone_index]+6*B+2] += (G-Rb)*S;
469 f[zofs[zone_index]+6*B+3] += Hb*S;
470 f[zofs[zone_index]+6*B+4] += Hnb*S;
471 f[zofs[zone_index]+6*B+5] += Hpb*S;
472
473 S = 0.25*ptri->d[1]*ptri->edge_len[1] + 0.25*ptri->d[0]*ptri->edge_len[0];
474 PetscScalar Hnc = (Rc_AUG-G)*(Egc+1.5*kb*Tpc) - 1.5*kb*Tnc*(Rc_SHR+Rc_DIR-G) - 1.5*kb*nc*(Tnc-Tc)/tao_enc;
475 PetscScalar Hpc = (Rc_AUG-G)*(Egc+1.5*kb*Tnc) - 1.5*kb*Tpc*(Rc_SHR+Rc_DIR-G) - 1.5*kb*pc*(Tpc-Tc)/tao_epc;
476 PetscScalar Hc = Rc_SHR*(Egc+1.5*kb*Tnc+1.5*kb*Tpc) + 1.5*kb*nc*(Tnc-Tc)/tao_enc + 1.5*kb*pc*(Tpc-Tc)/tao_epc;
477 f[zofs[zone_index]+6*C+1] += (G-Rc)*S;
478 f[zofs[zone_index]+6*C+2] += (G-Rc)*S;
479 f[zofs[zone_index]+6*C+3] += Hc*S;
480 f[zofs[zone_index]+6*C+4] += Hnc*S;
481 f[zofs[zone_index]+6*C+5] += Hpc*S;
482
483 }
484
485
486
487 void SMCZone::J3E_Tri_ddm(Tri *ptri,PetscScalar *x,Mat *jtmp, vector<int> & zofs)
488 {
489 int z = zone_index;
490 int A = ptri->node[0];
491 int B = ptri->node[1];
492 int C = ptri->node[2];
493 PetscScalar xa = pzone->danode[ptri->node[0]].x;
494 PetscScalar xb = pzone->danode[ptri->node[1]].x;
495 PetscScalar xc = pzone->danode[ptri->node[2]].x;
496 PetscScalar ya = pzone->danode[ptri->node[0]].y;
497 PetscScalar yb = pzone->danode[ptri->node[1]].y;
498 PetscScalar yc = pzone->danode[ptri->node[2]].y;
499 PetscScalar tri_area = ptri->area;
500 PetscScalar Sax = (xc-xb)/ptri->edge_len[0];
501 PetscScalar Say = (yc-yb)/ptri->edge_len[0];
502 PetscScalar Sbx = (xa-xc)/ptri->edge_len[1];
503 PetscScalar Sby = (ya-yc)/ptri->edge_len[1];
504 PetscScalar Scx = (xb-xa)/ptri->edge_len[2];
505 PetscScalar Scy = (yb-ya)/ptri->edge_len[2];
506 PetscScalar Wab = ptri->d[1]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
507 PetscScalar Wac = ptri->d[2]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[1])*sin(ptri->angle[1]));
508 PetscScalar Wbc = ptri->d[2]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[0])*sin(ptri->angle[0]));
509 PetscScalar Wba = ptri->d[0]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
510 PetscScalar Wca = ptri->d[0]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[1])*sin(ptri->angle[1]));
511 PetscScalar Wcb = ptri->d[1]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[0])*sin(ptri->angle[0]));
512 PetscScalar Ma = -cos(ptri->angle[0]);
513 PetscScalar Mb = -cos(ptri->angle[1]);
514 PetscScalar Mc = -cos(ptri->angle[2]);
515
516 PetscScalar e = mt->e;
517 PetscScalar kb = mt->kb;
518 PetscScalar S;
519 //the indepedent variable number
520 adtl::AutoDScalar::numdir=18;
521 //synchronize with material database
522 mt->set_ad_num(adtl::AutoDScalar::numdir);
523
524 Mat6 J1,J2,J3;
525 AutoDScalar sn,sp,F[6];
526 PetscInt index1[6],index2[6],index3[6];
527 index1[0] = zofs[z]+6*A+0;
528 index1[1] = zofs[z]+6*A+1;
529 index1[2] = zofs[z]+6*A+2;
530 index1[3] = zofs[z]+6*A+3;
531 index1[4] = zofs[z]+6*A+4;
532 index1[5] = zofs[z]+6*A+5;
533
534 index2[0] = zofs[z]+6*B+0;
535 index2[1] = zofs[z]+6*B+1;
536 index2[2] = zofs[z]+6*B+2;
537 index2[3] = zofs[z]+6*B+3;
538 index2[4] = zofs[z]+6*B+4;
539 index2[5] = zofs[z]+6*B+5;
540
541 index3[0] = zofs[z]+6*C+0;
542 index3[1] = zofs[z]+6*C+1;
543 index3[2] = zofs[z]+6*C+2;
544 index3[3] = zofs[z]+6*C+3;
545 index3[4] = zofs[z]+6*C+4;
546 index3[5] = zofs[z]+6*C+5;
547 //---------------------------------------------------------------------------
548 AutoDScalar Va = x[zofs[zone_index]+6*A+0]; //potential of node A
549 AutoDScalar na = x[zofs[zone_index]+6*A+1]; //electron density of node A
550 AutoDScalar pa = x[zofs[zone_index]+6*A+2]; //hole density of node A
551 AutoDScalar Ta = x[zofs[zone_index]+6*A+3]; //Lattice Temperature of node A
552 AutoDScalar naTna = x[zofs[zone_index]+6*A+4]; //na*Tna of node A
553 AutoDScalar paTpa = x[zofs[zone_index]+6*A+5]; //pa*Tpa of node A
554 Va.setADValue (0,1.0);
555 na.setADValue (1,1.0);
556 pa.setADValue (2,1.0);
557 Ta.setADValue (3,1.0);
558 naTna.setADValue(4,1.0);
559 paTpa.setADValue(5,1.0);
560 AutoDScalar Tna = naTna/na; // get Tna
561 AutoDScalar Tpa = paTpa/pa; // get Tpa
562 mt->mapping(&pzone->danode[A],&aux[A],0);
563 AutoDScalar nia = mt->band->nie(Ta);
564 AutoDScalar Ra = mt->band->Recomb(pa,na,Ta);
565 AutoDScalar Ra_SHR = mt->band->R_SHR(pa,na,Ta);
566 AutoDScalar Ra_AUG = mt->band->R_Auger(pa,na,Ta);
567 AutoDScalar Ra_DIR = mt->band->R_Direct(pa,na,Ta);
568 #ifdef _APPROXIMATE_
569 PetscScalar Nca = mt->band->Nc(fs[A].T);
570 PetscScalar Nva = mt->band->Nv(fs[A].T);
571 PetscScalar Ega = mt->band->Eg(fs[A].T);
572 AutoDScalar Eca = -(e*Va + aux[A].affinity + mt->band->EgNarrowToEc(fs[A].T) + kb*fs[A].T*(log(Nca)-1.5*log(fs[A].T)));
573 AutoDScalar Eva = -(e*Va + aux[A].affinity - mt->band->EgNarrowToEv(fs[A].T) - kb*fs[A].T*(log(Nva)-1.5*log(fs[A].T)) + Ega);
574 #else
575 AutoDScalar Nca = mt->band->Nc(Ta);
576 AutoDScalar Nva = mt->band->Nv(Ta);
577 AutoDScalar Ega = mt->band->Eg(Ta);
578 AutoDScalar Eca = -(e*Va + aux[A].affinity + mt->band->EgNarrowToEc(Ta) + kb*Ta*(log(Nca)-1.5*log(Ta)));
579 AutoDScalar Eva = -(e*Va + aux[A].affinity - mt->band->EgNarrowToEv(Ta) - kb*Ta*(log(Nva)-1.5*log(Ta)) + Ega);
580 #endif
581 AutoDScalar tao_ena = mt->band->ElecEnergyRelaxTime(Tna,Ta);
582 AutoDScalar tao_epa = mt->band->HoleEnergyRelaxTime(Tpa,Ta);
583
584
585 AutoDScalar Vb = x[zofs[zone_index]+6*B+0]; //potential of node B
586 AutoDScalar nb = x[zofs[zone_index]+6*B+1]; //electron density of node B
587 AutoDScalar pb = x[zofs[zone_index]+6*B+2]; //hole density of node B
588 AutoDScalar Tb = x[zofs[zone_index]+6*B+3]; //Temperature of node B
589 AutoDScalar nbTnb = x[zofs[zone_index]+6*B+4]; //nb*Tnb of node B
590 AutoDScalar pbTpb = x[zofs[zone_index]+6*B+5]; //pb*Tpb of node B
591 Vb.setADValue (6,1.0);
592 nb.setADValue (7,1.0);
593 pb.setADValue (8,1.0);
594 Tb.setADValue (9,1.0);
595 nbTnb.setADValue(10,1.0);
596 pbTpb.setADValue(11,1.0);
597 AutoDScalar Tnb = nbTnb/nb; // get Tnb
598 AutoDScalar Tpb = pbTpb/pb; // get Tpb
599 mt->mapping(&pzone->danode[B],&aux[B],0);
600 AutoDScalar nib = mt->band->nie(Tb);
601 AutoDScalar Rb = mt->band->Recomb(pb,nb,Tb);
602 AutoDScalar Rb_SHR = mt->band->R_SHR(pb,nb,Tb);
603 AutoDScalar Rb_AUG = mt->band->R_Auger(pb,nb,Tb);
604 AutoDScalar Rb_DIR = mt->band->R_Direct(pb,nb,Tb);
605 #ifdef _APPROXIMATE_
606 PetscScalar Ncb = mt->band->Nc(fs[B].T);
607 PetscScalar Nvb = mt->band->Nv(fs[B].T);
608 PetscScalar Egb = mt->band->Eg(fs[B].T);
609 AutoDScalar Ecb = -(e*Vb + aux[B].affinity + mt->band->EgNarrowToEc(fs[B].T)+ kb*fs[B].T*(log(Ncb)-1.5*log(fs[B].T)));
610 AutoDScalar Evb = -(e*Vb + aux[B].affinity - mt->band->EgNarrowToEv(fs[B].T)- kb*fs[B].T*(log(Nvb)-1.5*log(fs[B].T))+ Egb);
611 #else
612 AutoDScalar Ncb = mt->band->Nc(Tb);
613 AutoDScalar Nvb = mt->band->Nv(Tb);
614 AutoDScalar Egb = mt->band->Eg(Tb);
615 AutoDScalar Ecb = -(e*Vb + aux[B].affinity + mt->band->EgNarrowToEc(Tb)+ kb*Tb*(log(Ncb)-1.5*log(Tb)));
616 AutoDScalar Evb = -(e*Vb + aux[B].affinity - mt->band->EgNarrowToEv(Tb)- kb*Tb*(log(Nvb)-1.5*log(Tb))+ Egb);
617 #endif
618 AutoDScalar tao_enb = mt->band->ElecEnergyRelaxTime(Tnb,Tb);
619 AutoDScalar tao_epb = mt->band->HoleEnergyRelaxTime(Tpb,Tb);
620
621 AutoDScalar Vc = x[zofs[zone_index]+6*C+0]; //potential of node C
622 AutoDScalar nc = x[zofs[zone_index]+6*C+1]; //electron density of node C
623 AutoDScalar pc = x[zofs[zone_index]+6*C+2]; //hole density of node C
624 AutoDScalar Tc = x[zofs[zone_index]+6*C+3]; //Temperature of node C
625 AutoDScalar ncTnc = x[zofs[zone_index]+6*C+4]; //nc*Tnc of node C
626 AutoDScalar pcTpc = x[zofs[zone_index]+6*C+5]; //pc*Tpc of node C
627 Vc.setADValue (12,1.0);
628 nc.setADValue (13,1.0);
629 pc.setADValue (14,1.0);
630 Tc.setADValue (15,1.0);
631 ncTnc.setADValue(16,1.0);
632 pcTpc.setADValue(17,1.0);
633 AutoDScalar Tnc = ncTnc/nc; // get Tnb
634 AutoDScalar Tpc = pcTpc/pc; // get Tpb
635 mt->mapping(&pzone->danode[C],&aux[C],0);
636 AutoDScalar nic = mt->band->nie(Tc);
637 AutoDScalar Rc = mt->band->Recomb(pc,nc,Tc);
638 AutoDScalar Rc_SHR = mt->band->R_SHR(pc,nc,Tc);
639 AutoDScalar Rc_AUG = mt->band->R_Auger(pc,nc,Tc);
640 AutoDScalar Rc_DIR = mt->band->R_Direct(pc,nc,Tc);
641 #ifdef _APPROXIMATE_
642 PetscScalar Ncc = mt->band->Nc(fs[C].T);
643 PetscScalar Nvc = mt->band->Nv(fs[C].T);
644 PetscScalar Egc = mt->band->Eg(fs[C].T);
645 AutoDScalar Ecc = -(e*Vc + aux[C].affinity + mt->band->EgNarrowToEc(fs[C].T)+ kb*fs[C].T*(log(Ncc)-1.5*log(fs[C].T)));
646 AutoDScalar Evc = -(e*Vc + aux[C].affinity - mt->band->EgNarrowToEv(fs[C].T)- kb*fs[C].T*(log(Nvc)-1.5*log(fs[C].T))+ Egc);
647 #else
648 AutoDScalar Ncc = mt->band->Nc(Tc);
649 AutoDScalar Nvc = mt->band->Nv(Tc);
650 AutoDScalar Egc = mt->band->Eg(Tc);
651 AutoDScalar Ecc = -(e*Vc + aux[C].affinity + mt->band->EgNarrowToEc(Tc)+ kb*Tc*(log(Ncc)-1.5*log(Tc)));
652 AutoDScalar Evc = -(e*Vc + aux[C].affinity - mt->band->EgNarrowToEv(Tc)- kb*Tc*(log(Nvc)-1.5*log(Tc))+ Egc);
653 #endif
654 AutoDScalar tao_enc = mt->band->ElecEnergyRelaxTime(Tnc,Tc);
655 AutoDScalar tao_epc = mt->band->HoleEnergyRelaxTime(Tpc,Tc);
656
657
658 AutoDScalar Ex = -((yb-yc)*Va + (yc-ya)*Vb +(ya-yb)*Vc)/(2*tri_area);
659 AutoDScalar Ey = -((xc-xb)*Va + (xa-xc)*Vb +(xb-xa)*Vc)/(2*tri_area);
660 AutoDScalar Eg = mt->band->Eg((Ta+Tb+Tc)/3.0);
661 AutoDScalar Hn,Hp,Ga=0,Gb=0,Gc=0,G=0;
662 AutoDScalar IIn=0,IIp=0;
663 AutoDScalar mun,mup;
664 AutoDScalar Epnc=0,Etnc=0;
665 AutoDScalar Eppc=0,Etpc=0;
666 AutoDScalar Epna=0,Etna=0;
667 AutoDScalar Eppa=0,Etpa=0;
668 AutoDScalar Epnb=0,Etnb=0;
669 AutoDScalar Eppb=0,Etpb=0;
670 AutoDScalar Jna_norm,Jpa_norm;
671 AutoDScalar Jnb_norm,Jpb_norm;
672 AutoDScalar Jnc_norm,Jpc_norm;
673 AutoDScalar grad_P,kapa,grad_T,T_mid;
674
675 AutoDScalar Jnc = In(kb, e, -Eca/e, -Ecb/e, na, nb, Tna,Tnb, ptri->edge_len[2]);
676 AutoDScalar Jpc = Ip(kb, e, -Eva/e, -Evb/e, pa, pb, Tpa,Tpb, ptri->edge_len[2]);
677 AutoDScalar Jna = In(kb, e, -Ecb/e, -Ecc/e, nb, nc, Tnb,Tnc, ptri->edge_len[0]);
678 AutoDScalar Jpa = Ip(kb, e, -Evb/e, -Evc/e, pb, pc, Tpb,Tpc, ptri->edge_len[0]);
679 AutoDScalar Jnb = In(kb, e, -Ecc/e, -Eca/e, nc, na, Tnc,Tna, ptri->edge_len[1]);
680 AutoDScalar Jpb = Ip(kb, e, -Evc/e, -Eva/e, pc, pa, Tpc,Tpa, ptri->edge_len[1]);
681 PetscScalar Jn_scale = dmax(dmax(fabs(Jna.getValue()),fabs(Jnb.getValue())),
682 dmax(fabs(Jnc.getValue()),nia.getValue()*nia.getValue()));
683 PetscScalar Jp_scale = dmax(dmax(fabs(Jpa.getValue()),fabs(Jpb.getValue())),
684 dmax(fabs(Jpc.getValue()),nia.getValue()*nia.getValue()));
685 Jnc /= Jn_scale;
686 Jpc /= Jp_scale;
687 Jna /= Jn_scale;
688 Jpa /= Jp_scale;
689 Jnb /= Jn_scale;
690 Jpb /= Jp_scale;
691
692 //---------------------------------------------------------------------------
693 //process edge A-B
694 //---------------------------------------------------------------------------
695 if(HighFieldMobility)
696 {
697 if(EJModel || IIType==EdotJ || IIType==TempII)
698 {
699 AutoDScalar Jncx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jnc
700 + (Wca*Sax - Wca*Mb*Scx)*Jna
701 + (Wcb*Sbx - Wcb*Ma*Scx)*Jnb;
702 AutoDScalar Jncy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jnc
703 + (Wca*Say - Wca*Mb*Scy)*Jna
704 + (Wcb*Sby - Wcb*Ma*Scy)*Jnb;
705 AutoDScalar Jpcx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jpc
706 + (Wca*Sax - Wca*Mb*Scx)*Jpa
707 + (Wcb*Sbx - Wcb*Ma*Scx)*Jpb;
708 AutoDScalar Jpcy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jpc
709 + (Wca*Say - Wca*Mb*Scy)*Jpa
710 + (Wcb*Sby - Wcb*Ma*Scy)*Jpb;
711 Jnc_norm = sqrt(Jncx*Jncx + Jncy*Jncy + 1e-100);
712 Jpc_norm = sqrt(Jpcx*Jpcx + Jpcy*Jpcy + 1e-100);
713 Epnc = (Ex*Jncx + Ey*Jncy)/Jnc_norm;
714 Etnc = (Ex*Jncy - Ey*Jncx)/Jnc_norm;
715 Eppc = (Ex*Jpcx + Ey*Jpcy)/Jpc_norm;
716 Etpc = (Ex*Jpcy - Ey*Jpcx)/Jpc_norm;
717 }
718 else
719 {
720 Epnc = fabs(Eca-Ecb)/e/ptri->edge_len[2]; //parallel electrical field for electron
721 Eppc = fabs(Eva-Evb)/e/ptri->edge_len[2]; //parallel electrical field for hole
722 //transvers electrical field for electron and hole
723 Etnc = fabs(Eca + ptri->edge_len[1]*cos(ptri->angle[0])*(Ecb-Eca)/ptri->edge_len[2] - Ecc)
724 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
725 Etpc = fabs(Eva + ptri->edge_len[1]*cos(ptri->angle[0])*(Evb-Eva)/ptri->edge_len[2] - Evc)
726 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
727 }
728
729 mt->mapping(&pzone->danode[A],&aux[A],0);
730 AutoDScalar munCA = mt->mob->ElecMob(pa,na,Ta,adtl::fmax(0,Epnc),fabs(Etnc),1);
731 AutoDScalar mupCA = mt->mob->HoleMob(pa,na,Ta,adtl::fmax(0,Eppc),fabs(Etpc),1);
732
733 mt->mapping(&pzone->danode[B],&aux[B],0);
734 AutoDScalar munCB = mt->mob->ElecMob(pb,nb,Tb,adtl::fmax(0,Epnc),fabs(Etnc),1);
735 AutoDScalar mupCB = mt->mob->HoleMob(pb,nb,Tb,adtl::fmax(0,Eppc),fabs(Etpc),1);
736
737 mun = 0.5*(munCA+munCB);
738 mup = 0.5*(mupCA+mupCB);
739 }
740 else
741 {
742 mun = 0.5*(aux[A].mun+aux[B].mun);
743 mup = 0.5*(aux[A].mup+aux[B].mup);
744 }
745
746 T_mid = 0.5*(Ta+Tb);
747 grad_P = (Vb-Va)/ptri->edge_len[2];
748 kapa = mt->thermal->HeatConduction(T_mid);
749 grad_T = kapa*(Tb-Ta)/ptri->edge_len[2];
750 sn = Sn(kb, e, -Eca/e, -Ecb/e, na, nb, Tna,Tnb, ptri->edge_len[2]);
751 sp = Sp(kb, e, -Eva/e, -Evb/e, pa, pb, Tpa,Tpb, ptri->edge_len[2]);
752
753 //---------------------------------------------------------------------------
754 // even items
755 Hn = 0.5*(Va-Vb)*mun*Jnc*Jn_scale*ptri->d[2];
756 Hp = 0.5*(Va-Vb)*mup*Jpc*Jp_scale*ptri->d[2];
757 if(ImpactIonization)
758 {
759 if(IIType==EdotJ)
760 {
761 IIn = mt->gen->ElecGenRate(0.5*(Ta+Tb),adtl::fmax(0,Epnc),Eg);
762 IIp = mt->gen->HoleGenRate(0.5*(Ta+Tb),adtl::fmax(0,Eppc),Eg);
763 Gc = IIn*mun*Jnc_norm*Jn_scale + IIp*mup*Jpc_norm*Jp_scale;
764 }
765 if(IIType==TempII)
766 {
767 IIn = mt->gen->ElecGenRateEBM(0.5*(Tna+Tnb),0.5*(Ta+Tb),Eg);
768 IIp = mt->gen->HoleGenRateEBM(0.5*(Tpa+Tpb),0.5*(Ta+Tb),Eg);
769 Gc = IIn*mun*Jnc_norm*Jn_scale + IIp*mup*Jpc_norm*Jp_scale;
770 }
771 if(IIType==ESide )
772 {
773 IIn = mt->gen->ElecGenRate(T_mid,fabs(Va-Vb)/ptri->edge_len[2],Eg);
774 IIp = mt->gen->HoleGenRate(T_mid,fabs(Va-Vb)/ptri->edge_len[2],Eg);
775 Gc = IIn*mun*fabs(Jnc)*Jn_scale + IIp*mup*fabs(Jpc)*Jp_scale;
776 }
777
778 }
779 Set_Mat6_zero(J1);
780 Set_Mat6_zero(J2);
781 Set_Mat6_zero(J3);
782 for(int i=0;i<6;i++)
783 {
784 J1.m[1*6+i] = Gc.getADValue(0+i)*0.25*ptri->d[2]*ptri->edge_len[2];
785 J1.m[2*6+i] = Gc.getADValue(0+i)*0.25*ptri->d[2]*ptri->edge_len[2];
786 J1.m[4*6+i] = Hn.getADValue(0+i);
787 J1.m[5*6+i] = Hp.getADValue(0+i);
788
789 J2.m[1*6+i] = Gc.getADValue(6+i)*0.25*ptri->d[2]*ptri->edge_len[2];
790 J2.m[2*6+i] = Gc.getADValue(6+i)*0.25*ptri->d[2]*ptri->edge_len[2];
791 J2.m[4*6+i] = Hn.getADValue(6+i);
792 J2.m[5*6+i] = Hp.getADValue(6+i);
793
794 J3.m[1*6+i] = Gc.getADValue(12+i)*0.25*ptri->d[2]*ptri->edge_len[2];
795 J3.m[2*6+i] = Gc.getADValue(12+i)*0.25*ptri->d[2]*ptri->edge_len[2];
796 J3.m[4*6+i] = Hn.getADValue(12+i);
797 J3.m[5*6+i] = Hp.getADValue(12+i);
798 }
799 MatSetValues(*jtmp,6,index1,6,index1,J1.m,ADD_VALUES);
800 MatSetValues(*jtmp,6,index1,6,index2,J2.m,ADD_VALUES);
801 MatSetValues(*jtmp,6,index1,6,index3,J3.m,ADD_VALUES);
802 MatSetValues(*jtmp,6,index2,6,index1,J1.m,ADD_VALUES);
803 MatSetValues(*jtmp,6,index2,6,index2,J2.m,ADD_VALUES);
804 MatSetValues(*jtmp,6,index2,6,index3,J3.m,ADD_VALUES);
805 //---------------------------------------------------------------------------
806 // odd items
807 F[0] = grad_P*ptri->d[2];
808 F[1] = mun*Jnc*Jn_scale*ptri->d[2];
809 F[2] = - mup*Jpc*Jp_scale*ptri->d[2];
810 F[3] = grad_T*ptri->d[2];
811 F[4] = - mun*sn*ptri->d[2];
812 F[5] = - mup*sp*ptri->d[2];
813 for(int i=0;i<6;i++)
814 for(int j=0;j<6;j++)
815 {
816 J1.m[6*i+j] = F[i].getADValue(0+j);
817 J2.m[6*i+j] = F[i].getADValue(6+j);
818 J3.m[6*i+j] = F[i].getADValue(12+j);
819 }
820
821 MatSetValues(*jtmp,6,index1,6,index1,J1.m,ADD_VALUES);
822 MatSetValues(*jtmp,6,index1,6,index2,J2.m,ADD_VALUES);
823 MatSetValues(*jtmp,6,index1,6,index3,J3.m,ADD_VALUES);
824 MatSetValues(*jtmp,6,index2,6,index1,(-J1).m,ADD_VALUES);
825 MatSetValues(*jtmp,6,index2,6,index2,(-J2).m,ADD_VALUES);
826 MatSetValues(*jtmp,6,index2,6,index3,(-J3).m,ADD_VALUES);
827
828
829 //---------------------------------------------------------------------------
830 //process edge B-C
831 //---------------------------------------------------------------------------
832 if(HighFieldMobility)
833 {
834 if(EJModel || IIType==EdotJ || IIType==TempII)
835 {
836 AutoDScalar Jnax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jna
837 + (Wab*Sbx - Wab*Mc*Sax)*Jnb
838 + (Wac*Scx - Wac*Mb*Sax)*Jnc;
839 AutoDScalar Jnay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jna
840 + (Wab*Sby - Wab*Mc*Say)*Jnb
841 + (Wac*Scy - Wac*Mb*Say)*Jnc;
842 AutoDScalar Jpax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jpa
843 + (Wab*Sbx - Wab*Mc*Sax)*Jpb
844 + (Wac*Scx - Wac*Mb*Sax)*Jpc;
845 AutoDScalar Jpay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jpa
846 + (Wab*Sby - Wab*Mc*Say)*Jpb
847 + (Wac*Scy - Wac*Mb*Say)*Jpc;
848 Jna_norm = sqrt(Jnax*Jnax + Jnay*Jnay + 1e-100);
849 Jpa_norm = sqrt(Jpax*Jpax + Jpay*Jpay + 1e-100);
850 Epna = (Ex*Jnax + Ey*Jnay)/Jna_norm;
851 Etna = (Ex*Jnay - Ey*Jnax)/Jna_norm;
852 Eppa = (Ex*Jpax + Ey*Jpay)/Jpa_norm;
853 Etpa = (Ex*Jpay - Ey*Jpax)/Jpa_norm;
854 }
855 else
856 {
857 Epna = fabs(Ecb-Ecc)/e/ptri->edge_len[0]; //parallel electrical field for electron
858 Eppa = fabs(Evb-Evc)/e/ptri->edge_len[0]; //parallel electrical field for hole
859 //transvers electrical field for electron and hole
860 Etna = fabs(Ecb + ptri->edge_len[2]*cos(ptri->angle[1])*(Ecc-Ecb)/ptri->edge_len[0] - Eca)
861 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
862 Etpa = fabs(Evb + ptri->edge_len[2]*cos(ptri->angle[1])*(Evc-Evb)/ptri->edge_len[0] - Eva)
863 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
864 }
865 mt->mapping(&pzone->danode[B],&aux[B],0);
866 AutoDScalar munAB = mt->mob->ElecMob(pb,nb,Tb,adtl::fmax(0,Epna),fabs(Etna),1);
867 AutoDScalar mupAB = mt->mob->HoleMob(pb,nb,Tb,adtl::fmax(0,Eppa),fabs(Etpa),1);
868
869 mt->mapping(&pzone->danode[C],&aux[C],0);
870 AutoDScalar munAC = mt->mob->ElecMob(pc,nc,Tc,adtl::fmax(0,Epna),fabs(Etna),1);
871 AutoDScalar mupAC = mt->mob->HoleMob(pc,nc,Tc,adtl::fmax(0,Eppa),fabs(Etpa),1);
872
873 mun = 0.5*(munAB+munAC);
874 mup = 0.5*(mupAB+mupAC);
875 }
876 else
877 {
878 mun = 0.5*(aux[B].mun+aux[C].mun);
879 mup = 0.5*(aux[B].mup+aux[C].mup);
880 }
881 T_mid = 0.5*(Tb+Tc);
882 grad_P = (Vc-Vb)/ptri->edge_len[0];
883 kapa = mt->thermal->HeatConduction(T_mid);
884 grad_T = kapa*(Tc-Tb)/ptri->edge_len[0];
885
886 sn = Sn(kb, e, -Ecb/e, -Ecc/e, nb, nc, Tnb,Tnc, ptri->edge_len[0]);
887 sp = Sp(kb, e, -Evb/e, -Evc/e, pb, pc, Tpb,Tpc, ptri->edge_len[0]);
888 //---------------------------------------------------------------------------
889 // even item
890 Hn = 0.5*(Vb-Vc)*mun*Jna*Jn_scale*ptri->d[0];
891 Hp = 0.5*(Vb-Vc)*mup*Jpa*Jp_scale*ptri->d[0];
892 if(ImpactIonization)
893 {
894 if(IIType==EdotJ)
895 {
896 IIn = mt->gen->ElecGenRate(0.5*(Tb+Tc),adtl::fmax(0,Epna),Eg);
897 IIp = mt->gen->HoleGenRate(0.5*(Tb+Tc),adtl::fmax(0,Eppa),Eg);
898 Ga = IIn*mun*Jna_norm*Jn_scale + IIp*mup*Jpa_norm*Jp_scale;
899 }
900 if(IIType==TempII)
901 {
902 IIn = mt->gen->ElecGenRateEBM(0.5*(Tnb+Tnc),0.5*(Tb+Tc),Eg);
903 IIp = mt->gen->HoleGenRateEBM(0.5*(Tpb+Tpc),0.5*(Tb+Tc),Eg);
904 Ga = IIn*mun*Jna_norm*Jn_scale + IIp*mup*Jpa_norm*Jp_scale;
905 }
906 if(IIType==ESide )
907 {
908 IIn = mt->gen->ElecGenRate(T_mid,fabs(Vb-Vc)/ptri->edge_len[0],Eg);
909 IIp = mt->gen->HoleGenRate(T_mid,fabs(Vb-Vc)/ptri->edge_len[0],Eg);
910 Ga = IIn*mun*fabs(Jna)*Jn_scale + IIp*mup*fabs(Jpa)*Jp_scale;
911 }
912 }
913 Set_Mat6_zero(J1);
914 Set_Mat6_zero(J2);
915 Set_Mat6_zero(J3);
916 for(int i=0;i<6;i++)
917 {
918 J1.m[1*6+i] = Ga.getADValue(0+i)*0.25*ptri->d[0]*ptri->edge_len[0];
919 J1.m[2*6+i] = Ga.getADValue(0+i)*0.25*ptri->d[0]*ptri->edge_len[0];
920 J1.m[4*6+i] = Hn.getADValue(0+i);
921 J1.m[5*6+i] = Hp.getADValue(0+i);
922
923 J2.m[1*6+i] = Ga.getADValue(6+i)*0.25*ptri->d[0]*ptri->edge_len[0];
924 J2.m[2*6+i] = Ga.getADValue(6+i)*0.25*ptri->d[0]*ptri->edge_len[0];
925 J2.m[4*6+i] = Hn.getADValue(6+i);
926 J2.m[5*6+i] = Hp.getADValue(6+i);
927
928 J3.m[1*6+i] = Ga.getADValue(12+i)*0.25*ptri->d[0]*ptri->edge_len[0];
929 J3.m[2*6+i] = Ga.getADValue(12+i)*0.25*ptri->d[0]*ptri->edge_len[0];
930 J3.m[4*6+i] = Hn.getADValue(12+i);
931 J3.m[5*6+i] = Hp.getADValue(12+i);
932 }
933 MatSetValues(*jtmp,6,index2,6,index1,J1.m,ADD_VALUES);
934 MatSetValues(*jtmp,6,index2,6,index2,J2.m,ADD_VALUES);
935 MatSetValues(*jtmp,6,index2,6,index3,J3.m,ADD_VALUES);
936 MatSetValues(*jtmp,6,index3,6,index1,J1.m,ADD_VALUES);
937 MatSetValues(*jtmp,6,index3,6,index2,J2.m,ADD_VALUES);
938 MatSetValues(*jtmp,6,index3,6,index3,J3.m,ADD_VALUES);
939 //---------------------------------------------------------------------------
940 // odd item
941 F[0] = grad_P*ptri->d[0];
942 F[1] = mun*Jna*Jn_scale*ptri->d[0];
943 F[2] = - mup*Jpa*Jp_scale*ptri->d[0];
944 F[3] = grad_T*ptri->d[0];
945 F[4] = - mun*sn*ptri->d[0];
946 F[5] = - mup*sp*ptri->d[0];
947
948 for(int i=0;i<6;i++)
949 for(int j=0;j<6;j++)
950 {
951 J1.m[6*i+j] = F[i].getADValue(0+j);
952 J2.m[6*i+j] = F[i].getADValue(6+j);
953 J3.m[6*i+j] = F[i].getADValue(12+j);
954 }
955
956 MatSetValues(*jtmp,6,index2,6,index1,J1.m,ADD_VALUES);
957 MatSetValues(*jtmp,6,index2,6,index2,J2.m,ADD_VALUES);
958 MatSetValues(*jtmp,6,index2,6,index3,J3.m,ADD_VALUES);
959 MatSetValues(*jtmp,6,index3,6,index1,(-J1).m,ADD_VALUES);
960 MatSetValues(*jtmp,6,index3,6,index2,(-J2).m,ADD_VALUES);
961 MatSetValues(*jtmp,6,index3,6,index3,(-J3).m,ADD_VALUES);
962
963 //---------------------------------------------------------------------------
964 //process edge C-A
965 //---------------------------------------------------------------------------
966 if(HighFieldMobility)
967 {
968 if(EJModel || IIType==EdotJ || IIType==TempII)
969 {
970 AutoDScalar Jnbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jnb
971 + (Wbc*Scx - Wbc*Ma*Sbx)*Jnc
972 + (Wba*Sax - Wba*Mc*Sbx)*Jna;
973 AutoDScalar Jnby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jnb
974 + (Wbc*Scy - Wbc*Ma*Sby)*Jnc
975 + (Wba*Say - Wba*Mc*Sby)*Jna;
976 AutoDScalar Jpbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jpb
977 + (Wbc*Scx - Wbc*Ma*Sbx)*Jpc
978 + (Wba*Sax - Wba*Mc*Sbx)*Jpa;
979 AutoDScalar Jpby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jpb
980 + (Wbc*Scy - Wbc*Ma*Sby)*Jpc
981 + (Wba*Say - Wba*Mc*Sby)*Jpa;
982 Jnb_norm = sqrt(Jnbx*Jnbx + Jnby*Jnby + 1e-100);
983 Jpb_norm = sqrt(Jpbx*Jpbx + Jpby*Jpby + 1e-100);
984 Epnb = (Ex*Jnbx + Ey*Jnby)/Jnb_norm;
985 Etnb = (Ex*Jnby - Ey*Jnbx)/Jnb_norm;
986 Eppb = (Ex*Jpbx + Ey*Jpby)/Jpb_norm;
987 Etpb = (Ex*Jpby - Ey*Jpbx)/Jpb_norm;
988 }
989 else
990 {
991 Epnb = fabs(Ecc-Eca)/e/ptri->edge_len[1]; //parallel electrical field for electron
992 Eppb = fabs(Evc-Eva)/e/ptri->edge_len[1]; //parallel electrical field for hole
993 //transvers electrical field for electron and hole
994 Etnb = fabs(Ecc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eca-Ecc)/ptri->edge_len[1] - Ecb)
995 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
996 Etpb = fabs(Evc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eva-Evc)/ptri->edge_len[1] - Evb)
997 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
998 }
999 mt->mapping(&pzone->danode[C],&aux[C],0);
1000 AutoDScalar munBC = mt->mob->ElecMob(pc,nc,Tc,adtl::fmax(0,Epnb),fabs(Etnb),1);
1001 AutoDScalar mupBC = mt->mob->HoleMob(pc,nc,Tc,adtl::fmax(0,Eppb),fabs(Etpb),1);
1002
1003 mt->mapping(&pzone->danode[A],&aux[A],0);
1004 AutoDScalar munBA = mt->mob->ElecMob(pa,na,Ta,adtl::fmax(0,Epnb),fabs(Etnb),1);
1005 AutoDScalar mupBA = mt->mob->HoleMob(pa,na,Ta,adtl::fmax(0,Eppb),fabs(Etpb),1);
1006
1007 mun = 0.5*(munBC+munBA);
1008 mup = 0.5*(mupBC+mupBA);
1009 }
1010 else
1011 {
1012 mun = 0.5*(aux[C].mun+aux[A].mun);
1013 mup = 0.5*(aux[C].mup+aux[A].mup);
1014 }
1015 T_mid = 0.5*(Tc+Ta);
1016 grad_P = (Va-Vc)/ptri->edge_len[1];
1017 kapa = mt->thermal->HeatConduction(T_mid);
1018 grad_T = kapa*(Ta-Tc)/ptri->edge_len[1];
1019
1020 sn = Sn(kb, e, -Ecc/e, -Eca/e, nc, na, Tnc,Tna, ptri->edge_len[1]);
1021 sp = Sp(kb, e, -Evc/e, -Eva/e, pc, pa, Tpc,Tpa, ptri->edge_len[1]);
1022 //---------------------------------------------------------------------------
1023 // even item
1024 Hn = 0.5*(Vc-Va)*mun*Jnb*Jn_scale*ptri->d[1];
1025 Hp = 0.5*(Vc-Va)*mup*Jpb*Jp_scale*ptri->d[1];
1026 if(ImpactIonization)
1027 {
1028 if(IIType==EdotJ)
1029 {
1030 IIn = mt->gen->ElecGenRate(0.5*(Tc+Ta),adtl::fmax(0,Epnb),Eg);
1031 IIp = mt->gen->HoleGenRate(0.5*(Tc+Ta),adtl::fmax(0,Eppb),Eg);
1032 Gb = IIn*mun*Jnb_norm*Jn_scale + IIp*mup*Jpb_norm*Jp_scale;
1033 }
1034 if(IIType==TempII)
1035 {
1036 IIn = mt->gen->ElecGenRateEBM(0.5*(Tnc+Tna),0.5*(Tc+Ta),Eg);
1037 IIp = mt->gen->HoleGenRateEBM(0.5*(Tpc+Tpa),0.5*(Tc+Ta),Eg);
1038 Gb = IIn*mun*Jnb_norm*Jn_scale + IIp*mup*Jpb_norm*Jp_scale;
1039 }
1040 if(IIType==ESide )
1041 {
1042 IIn = mt->gen->ElecGenRate(T_mid,fabs(Vc-Va)/ptri->edge_len[1],Eg);
1043 IIp = mt->gen->HoleGenRate(T_mid,fabs(Vc-Va)/ptri->edge_len[1],Eg);
1044 Gb = IIn*mun*fabs(Jnb)*Jn_scale + IIp*mup*fabs(Jpb)*Jp_scale;
1045 }
1046 }
1047 Set_Mat6_zero(J1);
1048 Set_Mat6_zero(J2);
1049 Set_Mat6_zero(J3);
1050 for(int i=0;i<6;i++)
1051 {
1052 J1.m[1*6+i] = Gb.getADValue(0+i)*0.25*ptri->d[1]*ptri->edge_len[1];
1053 J1.m[2*6+i] = Gb.getADValue(0+i)*0.25*ptri->d[1]*ptri->edge_len[1];
1054 J1.m[4*6+i] = Hn.getADValue(0+i);
1055 J1.m[5*6+i] = Hp.getADValue(0+i);
1056
1057 J2.m[1*6+i] = Gb.getADValue(6+i)*0.25*ptri->d[1]*ptri->edge_len[1];
1058 J2.m[2*6+i] = Gb.getADValue(6+i)*0.25*ptri->d[1]*ptri->edge_len[1];
1059 J2.m[4*6+i] = Hn.getADValue(6+i);
1060 J2.m[5*6+i] = Hp.getADValue(6+i);
1061
1062 J3.m[1*6+i] = Gb.getADValue(12+i)*0.25*ptri->d[1]*ptri->edge_len[1];
1063 J3.m[2*6+i] = Gb.getADValue(12+i)*0.25*ptri->d[1]*ptri->edge_len[1];
1064 J3.m[4*6+i] = Hn.getADValue(12+i);
1065 J3.m[5*6+i] = Hp.getADValue(12+i);
1066 }
1067 MatSetValues(*jtmp,6,index3,6,index1,J1.m,ADD_VALUES);
1068 MatSetValues(*jtmp,6,index3,6,index2,J2.m,ADD_VALUES);
1069 MatSetValues(*jtmp,6,index3,6,index3,J3.m,ADD_VALUES);
1070 MatSetValues(*jtmp,6,index1,6,index1,J1.m,ADD_VALUES);
1071 MatSetValues(*jtmp,6,index1,6,index2,J2.m,ADD_VALUES);
1072 MatSetValues(*jtmp,6,index1,6,index3,J3.m,ADD_VALUES);
1073 //---------------------------------------------------------------------------
1074 // odd item
1075 F[0] = grad_P*ptri->d[1];
1076 F[1] = mun*Jnb*Jn_scale*ptri->d[1];
1077 F[2] = - mup*Jpb*Jp_scale*ptri->d[1];
1078 F[3] = grad_T*ptri->d[1];
1079 F[4] = - mun*sn*ptri->d[1];
1080 F[5] = - mup*sp*ptri->d[1];
1081 //---------------------------------------------------------------------------
1082 for(int i=0;i<6;i++)
1083 for(int j=0;j<6;j++)
1084 {
1085 J1.m[6*i+j] = F[i].getADValue(0+j);
1086 J2.m[6*i+j] = F[i].getADValue(6+j);
1087 J3.m[6*i+j] = F[i].getADValue(12+j);
1088 }
1089
1090 MatSetValues(*jtmp,6,index3,6,index1,J1.m,ADD_VALUES);
1091 MatSetValues(*jtmp,6,index3,6,index2,J2.m,ADD_VALUES);
1092 MatSetValues(*jtmp,6,index3,6,index3,J3.m,ADD_VALUES);
1093 MatSetValues(*jtmp,6,index1,6,index1,(-J1).m,ADD_VALUES);
1094 MatSetValues(*jtmp,6,index1,6,index2,(-J2).m,ADD_VALUES);
1095 MatSetValues(*jtmp,6,index1,6,index3,(-J3).m,ADD_VALUES);
1096
1097
1098 //---------------------------------------------------------------------------
1099 // G-R item and heat source for each node
1100 //---------------------------------------------------------------------------
1101 G=(Ga+Gb+Gc)/3.0;
1102 S = 0.25*ptri->d[2]*ptri->edge_len[2] + 0.25*ptri->d[1]*ptri->edge_len[1];
1103 AutoDScalar Hna = (Ra_AUG-G)*(Ega+1.5*kb*Tpa) - 1.5*kb*Tna*(Ra_SHR+Ra_DIR-G) - 1.5*kb*na*(Tna-Ta)/tao_ena;
1104 AutoDScalar Hpa = (Ra_AUG-G)*(Ega+1.5*kb*Tna) - 1.5*kb*Tpa*(Ra_SHR+Ra_DIR-G) - 1.5*kb*pa*(Tpa-Ta)/tao_epa;
1105 AutoDScalar Ha = Ra_SHR*(Ega+1.5*kb*Tna+1.5*kb*Tpa) + 1.5*kb*na*(Tna-Ta)/tao_ena + 1.5*kb*pa*(Tpa-Ta)/tao_epa;
1106 F[0] = 0;
1107 F[1] = (G-Ra)*S;
1108 F[2] = (G-Ra)*S;
1109 F[3] = Ha*S;
1110 F[4] = Hna*S;
1111 F[5] = Hpa*S;
1112 for(int i=0;i<6;i++)
1113 for(int j=0;j<6;j++)
1114 {
1115 J1.m[6*i+j] = F[i].getADValue(0+j);
1116 }
1117 MatSetValues(*jtmp,6,index1,6,index1,J1.m,ADD_VALUES);
1118
1119 S = 0.25*ptri->d[0]*ptri->edge_len[0] + 0.25*ptri->d[2]*ptri->edge_len[2];
1120 AutoDScalar Hnb = (Rb_AUG-G)*(Egb+1.5*kb*Tpb) - 1.5*kb*Tnb*(Rb_SHR+Rb_DIR-G) - 1.5*kb*nb*(Tnb-Tb)/tao_enb;
1121 AutoDScalar Hpb = (Rb_AUG-G)*(Egb+1.5*kb*Tnb) - 1.5*kb*Tpb*(Rb_SHR+Rb_DIR-G) - 1.5*kb*pb*(Tpb-Tb)/tao_epb;
1122 AutoDScalar Hb = Rb_SHR*(Egb+1.5*kb*Tnb+1.5*kb*Tpb) + 1.5*kb*nb*(Tnb-Tb)/tao_enb + 1.5*kb*pb*(Tpb-Tb)/tao_epb;
1123 F[0] = 0;
1124 F[1] = (G-Rb)*S;
1125 F[2] = (G-Rb)*S;
1126 F[3] = Hb*S;
1127 F[4] = Hnb*S;
1128 F[5] = Hpb*S;
1129 for(int i=0;i<6;i++)
1130 for(int j=0;j<6;j++)
1131 {
1132 J2.m[6*i+j] = F[i].getADValue(6+j);
1133 }
1134 MatSetValues(*jtmp,6,index2,6,index2,J2.m,ADD_VALUES);
1135
1136 S = 0.25*ptri->d[1]*ptri->edge_len[1] + 0.25*ptri->d[0]*ptri->edge_len[0];
1137 AutoDScalar Hnc = (Rc_AUG-G)*(Egc+1.5*kb*Tpc) - 1.5*kb*Tnc*(Rc_SHR+Rc_DIR-G) - 1.5*kb*nc*(Tnc-Tc)/tao_enc;
1138 AutoDScalar Hpc = (Rc_AUG-G)*(Egc+1.5*kb*Tnc) - 1.5*kb*Tpc*(Rc_SHR+Rc_DIR-G) - 1.5*kb*pc*(Tpc-Tc)/tao_epc;
1139 AutoDScalar Hc = Rc_SHR*(Egc+1.5*kb*Tnc+1.5*kb*Tpc) + 1.5*kb*nc*(Tnc-Tc)/tao_enc + 1.5*kb*pc*(Tpc-Tc)/tao_epc;
1140 F[0] = 0;
1141 F[1] = (G-Rc)*S;
1142 F[2] = (G-Rc)*S;
1143 F[3] = Hc*S;
1144 F[4] = Hnc*S;
1145 F[5] = Hpc*S;
1146 for(int i=0;i<6;i++)
1147 for(int j=0;j<6;j++)
1148 {
1149 J3.m[6*i+j] = F[i].getADValue(12+j);
1150 }
1151 MatSetValues(*jtmp,6,index3,6,index3,J3.m,ADD_VALUES);
1152
1153 }
1154
1155
1156 //-----------------------------------------------------------------------------
1157 // boundaries
1158 //-----------------------------------------------------------------------------
1159
1160 void SMCZone::F3E_ddm_inner(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs)
1161 {
1162 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1163 PetscScalar kb = mt->kb;
1164 PetscScalar e = mt->e;
1165 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1166 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1167 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1168 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1169 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1170 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1171
1172 f[zofs[zone_index]+6*i+0] = f[zofs[zone_index]+6*i+0]/pcell->area + e/aux[i].eps*((pi-aux[i].Na)+(aux[i].Nd-ni));
1173 f[zofs[zone_index]+6*i+1] = f[zofs[zone_index]+6*i+1]/pcell->area;
1174 f[zofs[zone_index]+6*i+2] = f[zofs[zone_index]+6*i+2]/pcell->area;
1175 f[zofs[zone_index]+6*i+3] = f[zofs[zone_index]+6*i+3]/pcell->area;
1176 f[zofs[zone_index]+6*i+4] = f[zofs[zone_index]+6*i+4]/pcell->area;
1177 f[zofs[zone_index]+6*i+5] = f[zofs[zone_index]+6*i+5]/pcell->area;
1178
1179 if(ODE_F.TimeDependent==true)
1180 {
1181 //second order
1182 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1183 {
1184 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1185 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1186 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1187 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1188 PetscScalar Wnt = (2-r)/(1-r)*1.5*ni*kb*Tn-1.0/(r*(1-r))*1.5*fs[i].n*kb*fs[i].Tn+(1-r)/r*1.5*fs[i].n_last*kb*fs[i].Tn_last;
1189 PetscScalar Wpt = (2-r)/(1-r)*1.5*pi*kb*Tp-1.0/(r*(1-r))*1.5*fs[i].p*kb*fs[i].Tp+(1-r)/r*1.5*fs[i].p_last*kb*fs[i].Tp_last;
1190 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1191 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1192 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1193 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1194 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1195 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1196 }
1197 else //first order
1198 {
1199 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1200 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1201 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1202 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1203 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1204 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1205 }
1206 }
1207
1208 }
1209
1210
1211
1212 void SMCZone::F3E_ddm_neumannbc(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs,DABC &bc)
1213 {
1214 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1215 NeumannBC *pbc = dynamic_cast <NeumannBC * >(bc.Get_pointer(pcell->bc_index-1));
1216 PetscScalar kb = mt->kb;
1217 PetscScalar e = mt->e;
1218 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1219 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1220 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1221 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1222 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1223 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1224
1225 PetscScalar grad_T=0;
1226 for(int j=0;j<pcell->nb_num;j++)
1227 {
1228 int nb = pcell->nb_array[j];
1229 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1230 if(j==0||j==pcell->nb_num-1)
1231 {
1232 PetscScalar h = pbc->heat_transfer;
1233 PetscScalar r = h*pbc->T_external;
1234 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj));
1235 }
1236 }
1237 f[zofs[zone_index]+6*i+0] = f[zofs[zone_index]+6*i+0]/pcell->area + e/aux[i].eps*((pi-aux[i].Na)+(aux[i].Nd-ni));
1238 f[zofs[zone_index]+6*i+1] = f[zofs[zone_index]+6*i+1]/pcell->area;
1239 f[zofs[zone_index]+6*i+2] = f[zofs[zone_index]+6*i+2]/pcell->area;
1240 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T)/pcell->area;
1241 f[zofs[zone_index]+6*i+4] = f[zofs[zone_index]+6*i+4]/pcell->area;
1242 f[zofs[zone_index]+6*i+5] = f[zofs[zone_index]+6*i+5]/pcell->area;
1243
1244 if(ODE_F.TimeDependent==true)
1245 {
1246 //second order
1247 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1248 {
1249 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1250 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1251 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1252 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1253 PetscScalar Wnt = (2-r)/(1-r)*1.5*ni*kb*Tn-1.0/(r*(1-r))*1.5*fs[i].n*kb*fs[i].Tn+(1-r)/r*1.5*fs[i].n_last*kb*fs[i].Tn_last;
1254 PetscScalar Wpt = (2-r)/(1-r)*1.5*pi*kb*Tp-1.0/(r*(1-r))*1.5*fs[i].p*kb*fs[i].Tp+(1-r)/r*1.5*fs[i].p_last*kb*fs[i].Tp_last;
1255 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1256 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1257 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1258 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1259 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1260 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1261 }
1262 else //first order
1263 {
1264 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1265 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1266 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1267 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1268 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1269 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1270 }
1271 }
1272 }
1273
1274
1275 void SMCZone::F3E_ddm_ombc_segment(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
1276 {
1277 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1278 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
1279 int z = zone_index;
1280 int equ_num = 6;
1281 int size = pzone->davcell.size();
1282 PetscScalar kb = mt->kb;
1283 PetscScalar Na = aux[i].Na;
1284 PetscScalar Nd = aux[i].Nd;
1285 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1286 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1287 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1288 PetscScalar Ti = x[zofs[z]+6*i+3];
1289 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1290 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1291 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
1292 PetscScalar nie = mt->band->nie(fs[i].T);
1293 PetscScalar Nc = mt->band->Nc(fs[i].T);
1294 PetscScalar Nv = mt->band->Nv(fs[i].T);
1295 PetscScalar electron_density,hole_density;
1296 int om_equ;
1297 for(int j=0;j<electrode.size();j++)
1298 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
1299
1300 f[zofs[z]+6*i+0] = Vi - kb*fs[i].T/mt->e*asinh((Nd-Na)/(2*nie)) + mt->kb*fs[i].T/2/mt->e*log(Nc/Nv)
1301 + mt->band->Eg(fs[i].T)/2/mt->e
1302 + aux[i].affinity -x[zofs[z]+equ_num*size+om_equ];
1303
1304 if(Na>Nd) //p-type
1305 {
1306 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1307 electron_density = nie*nie/hole_density;
1308 }
1309 else //n-type
1310 {
1311 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1312 hole_density = nie*nie/electron_density;
1313 }
1314 f[zofs[z]+6*i+1] = x[zofs[z]+6*i+1] - electron_density; //electron density
1315 f[zofs[z]+6*i+2] = x[zofs[z]+6*i+2] - hole_density; //hole density
1316
1317 PetscScalar grad_T=0;
1318 for(int j=0;j<pcell->nb_num;j++)
1319 {
1320 int nb = pcell->nb_array[j];
1321 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1322 PetscScalar T_mid = 0.5*(Tj+Ti);
1323 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
1324 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
1325 grad_T += kapa*pcell->elen[j]*dTdx_mid/pcell->area;
1326 if(j==0||j==pcell->nb_num-1)
1327 {
1328 PetscScalar h = pbc->heat_transfer;
1329 PetscScalar r = h*pbc->T_external;
1330 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj))/pcell->area;
1331 }
1332 }
1333 f[zofs[zone_index]+6*i+3] = grad_T;
1334 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1335 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1336
1337 if(ODE_F.TimeDependent==true)
1338 {
1339 //second order
1340 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1341 {
1342 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1343 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1344 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1345 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1346 }
1347 else //first order
1348 {
1349 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1350 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1351 }
1352 }
1353 }
1354
1355
1356 void SMCZone::F3E_ddm_ombc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1357 ElZone *pz, int n)
1358 {
1359 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1360 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
1361 int z = zone_index;
1362 int equ_num = 6;
1363 int size = pzone->davcell.size();
1364 PetscScalar kb = mt->kb;
1365 PetscScalar Na = aux[i].Na;
1366 PetscScalar Nd = aux[i].Nd;
1367 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1368 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1369 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1370 PetscScalar Ti = x[zofs[z]+6*i+3];
1371 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1372 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1373 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
1374 PetscScalar nie = mt->band->nie(fs[i].T);
1375 PetscScalar Nc = mt->band->Nc(fs[i].T);
1376 PetscScalar Nv = mt->band->Nv(fs[i].T);
1377 PetscScalar electron_density,hole_density;
1378 int om_equ;
1379 for(int j=0;j<electrode.size();j++)
1380 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
1381
1382 f[zofs[z]+6*i+0] = Vi - kb*fs[i].T/mt->e*asinh((Nd-Na)/(2*nie)) + mt->kb*fs[i].T/2/mt->e*log(Nc/Nv)
1383 + mt->band->Eg(fs[i].T)/2/mt->e
1384 + aux[i].affinity -x[zofs[z]+equ_num*size+om_equ];
1385
1386 if(Na>Nd) //p-type
1387 {
1388 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1389 electron_density = nie*nie/hole_density;
1390 }
1391 else //n-type
1392 {
1393 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1394 hole_density = nie*nie/electron_density;
1395 }
1396 f[zofs[z]+6*i+1] = x[zofs[z]+6*i+1] - electron_density; //electron density
1397 f[zofs[z]+6*i+2] = x[zofs[z]+6*i+2] - hole_density; //hole density
1398
1399 PetscScalar grad_T=0;
1400 for(int j=0;j<pcell->nb_num;j++)
1401 {
1402 int nb = pcell->nb_array[j];
1403 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1404 PetscScalar T_mid = 0.5*(Tj+Ti);
1405 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
1406 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
1407 grad_T += kapa*pcell->elen[j]*dTdx_mid;
1408 }
1409 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1410 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
1411 for(int j=0;j<ncell->nb_num;j++)
1412 {
1413 int nb = ncell->nb_array[j];
1414 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
1415 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
1416 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
1417 }
1418
1419 f[zofs[zone_index]+6*i+3] = grad_T;
1420 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1421 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1422
1423
1424 if(ODE_F.TimeDependent==true)
1425 {
1426 //second order
1427 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1428 {
1429 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1430 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1431 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1432 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1433 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
1434 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
1435 }
1436 else //first order
1437 {
1438 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1439 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1440 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
1441 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
1442 }
1443 }
1444
1445 }
1446
1447
1448
1449 void SMCZone::F3E_ddm_stkbc_segment(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc)
1450 {
1451 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1452 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
1453 int z = zone_index;
1454 int equ_num = 6;
1455 int size = pzone->davcell.size();
1456 PetscScalar kb = mt->kb;
1457 PetscScalar e = mt->e;
1458 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1459 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1460 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1461 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
1462 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1463 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1464
1465 int stk_equ;
1466 for(int j=0;j<electrode.size();j++)
1467 if(electrode[j]==pcell->bc_index-1)
1468 {
1469 stk_equ=j;
1470 break;
1471 }
1472 //Schotty Barrier Lowerring
1473 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
1474 //potential
1475 f[zofs[z]+6*i+0] = x[zofs[z]+6*i+0] + pbc->WorkFunction - deltaVB - x[zofs[z]+equ_num*size+stk_equ];
1476
1477 PetscScalar Fn=0,Fp=0,grad_T=0;
1478 for(int j=0;j<pcell->nb_num;j++)
1479 {
1480 int nb = pcell->nb_array[j];
1481 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1482 if(j==0||j==pcell->nb_num-1)
1483 {
1484 Fn += mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[i].affinity - deltaVB)*0.5*pcell->ilen[j];
1485 Fp += mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[i].affinity + deltaVB)*0.5*pcell->ilen[j];
1486 PetscScalar h = pbc->heat_transfer;
1487 PetscScalar r = h*pbc->T_external;
1488 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj));
1489 }
1490 }
1491
1492 f[zofs[zone_index]+6*i+1] = (f[zofs[zone_index]+6*i+1]+Fn)/pcell->area;
1493 f[zofs[zone_index]+6*i+2] = (f[zofs[zone_index]+6*i+2]-Fp)/pcell->area;
1494 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T)/pcell->area;
1495 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1496 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1497
1498 if(ODE_F.TimeDependent==true)
1499 {
1500 //second order
1501 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1502 {
1503 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1504 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1505 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1506 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1507 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1508 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1509 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1510 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1511 }
1512 else //first order
1513 {
1514 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1515 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1516 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1517 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1518 }
1519 }
1520 }
1521
1522
1523 void SMCZone::F3E_ddm_stkbc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc,
1524 ElZone *pz, int n)
1525 {
1526 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1527 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
1528 int z = zone_index;
1529 int equ_num = 6;
1530 int size = pzone->davcell.size();
1531 PetscScalar kb = mt->kb;
1532 PetscScalar e = mt->e;
1533 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1534 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1535 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1536 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
1537 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1538 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1539
1540 int stk_equ;
1541 for(int j=0;j<electrode.size();j++)
1542 if(electrode[j]==pcell->bc_index-1)
1543 {
1544 stk_equ=j;
1545 break;
1546 }
1547 //Schotty Barrier Lowerring
1548 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
1549 //potential
1550 f[zofs[z]+6*i+0] = x[zofs[z]+6*i+0] + pbc->WorkFunction - deltaVB - x[zofs[z]+equ_num*size+stk_equ];
1551
1552 PetscScalar Fn=0,Fp=0,grad_T=0;
1553 for(int j=0;j<pcell->nb_num;j++)
1554 {
1555 int nb = pcell->nb_array[j];
1556 if(j==0||j==pcell->nb_num-1)
1557 {
1558 Fn += mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[i].affinity - deltaVB)*0.5*pcell->ilen[j];
1559 Fp += mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[i].affinity + deltaVB)*0.5*pcell->ilen[j];
1560 }
1561 }
1562
1563 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1564 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
1565 for(int j=0;j<ncell->nb_num;j++)
1566 {
1567 int nb = ncell->nb_array[j];
1568 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
1569 PetscScalar kapa = pz->mt->thermal->HeatConduction(fs[i].T);//(0.5*(Ti+Tj_n));
1570 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
1571 }
1572
1573 f[zofs[zone_index]+6*i+1] = (f[zofs[zone_index]+6*i+1]+Fn)/pcell->area;
1574 f[zofs[zone_index]+6*i+2] = (f[zofs[zone_index]+6*i+2]-Fp)/pcell->area;
1575 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T);
1576 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1577 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1578
1579 if(ODE_F.TimeDependent==true)
1580 {
1581 //second order
1582 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1583 {
1584 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1585 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1586 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1587 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1588 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1589 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1590 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1591 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1592 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
1593 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
1594 }
1595 else //first order
1596 {
1597 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1598 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1599 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1600 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1601 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
1602 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
1603 }
1604 }
1605
1606 }
1607
1608
1609 void SMCZone::F3E_ddm_insulator_gate(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
1610 {
1611 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1612 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
1613 int equ_num = 6;
1614 int size = pzone->davcell.size();
1615 PetscScalar kb = mt->kb;
1616 PetscScalar e = mt->e;
1617 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1618 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1619 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1620 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1621 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1622 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1623 PetscScalar grad_P = 0, grad_T = 0;
1624 int ins_equ;
1625 for(int j=0;j<electrode.size();j++)
1626 if(electrode[j]==pcell->bc_index-1)
1627 {ins_equ=j;break;}
1628 for(int j=0;j<pcell->nb_num;j++)
1629 {
1630 int nb = pcell->nb_array[j];
1631 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
1632 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1633
1634 if(j==0||j==pcell->nb_num-1)
1635 {
1636 PetscScalar vgate = x[zofs[zone_index]+equ_num*size+ins_equ] - pbc->WorkFunction;
1637 PetscScalar q = e*pbc->QF; //sigma is the surface change density
1638 PetscScalar Thick = pbc->Thick;
1639 PetscScalar eps_ox = mt->eps0*pbc->eps;
1640 PetscScalar eps = aux[i].eps;
1641 PetscScalar r=q/eps + eps_ox/eps/Thick*vgate;
1642 PetscScalar s=eps_ox/eps/Thick;
1643 grad_P += 0.5*pcell->ilen[j]*(r-0.25*s*(3*Vi+Vj));
1644
1645 PetscScalar h = pbc->heat_transfer;
1646 grad_T += 0.5*pcell->ilen[j]*(h*pbc->T_external-0.25*h*(3*Ti+Tj));
1647 }
1648 }
1649
1650 f[zofs[zone_index]+6*i+0] = (f[zofs[zone_index]+6*i+0]+grad_P)/pcell->area
1651 + e/aux[i].eps*((pi-aux[i].Na)+(aux[i].Nd-ni));
1652 f[zofs[zone_index]+6*i+1] = f[zofs[zone_index]+6*i+1]/pcell->area;
1653 f[zofs[zone_index]+6*i+2] = f[zofs[zone_index]+6*i+2]/pcell->area;
1654 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T)/pcell->area;
1655 f[zofs[zone_index]+6*i+4] = f[zofs[zone_index]+6*i+4]/pcell->area;
1656 f[zofs[zone_index]+6*i+5] = f[zofs[zone_index]+6*i+5]/pcell->area;
1657
1658 if(ODE_F.TimeDependent==true)
1659 {
1660 //second order
1661 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1662 {
1663 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1664 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1665 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1666 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1667 PetscScalar Wnt = (2-r)/(1-r)*1.5*ni*kb*Tn-1.0/(r*(1-r))*1.5*fs[i].n*kb*fs[i].Tn+(1-r)/r*1.5*fs[i].n_last*kb*fs[i].Tn_last;
1668 PetscScalar Wpt = (2-r)/(1-r)*1.5*pi*kb*Tp-1.0/(r*(1-r))*1.5*fs[i].p*kb*fs[i].Tp+(1-r)/r*1.5*fs[i].p_last*kb*fs[i].Tp_last;
1669 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1670 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1671 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1672 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1673 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1674 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1675 }
1676 else //first order
1677 {
1678 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1679 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1680 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1681 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1682 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1683 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1684 }
1685 }
1686 }
1687
1688
1689
1690 void SMCZone::F3E_ddm_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1691 ISZone *pz, int n)
1692 {
1693 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1694 InsulatorInterfaceBC *pbc = dynamic_cast<InsulatorInterfaceBC * >(bc.Get_pointer(pcell->bc_index-1));
1695 PetscScalar kb = mt->kb;
1696 PetscScalar e = mt->e;
1697 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1698 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1699 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1700 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1701 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1702 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1703
1704 PetscScalar grad_P = 0,grad_T=0;
1705 //semiconductor-insulator interface
1706 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1707 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
1708 for(int j=0;j<ncell->nb_num;j++)
1709 {
1710 int nb = ncell->nb_array[j];
1711 PetscScalar Vj_n = x[zofs[pz->pzone->zone_index]+2*nb+0]; //potential of nb node
1712 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
1713 grad_P += pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j]*(Vj_n-Vi);
1714 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
1715 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
1716 }
1717 PetscScalar L = (pcell->ilen[0]+pcell->ilen[pcell->nb_num-1])/2.0;
1718
1719
1720 f[zofs[zone_index]+6*i+0] = (aux[i].eps*f[zofs[zone_index]+6*i+0]+grad_P)
1721 + e*((pi-aux[i].Na)+(aux[i].Nd-ni))*pcell->area + pbc->QF*L;
1722 f[zofs[zone_index]+6*i+1] = f[zofs[zone_index]+6*i+1]/pcell->area;
1723 f[zofs[zone_index]+6*i+2] = f[zofs[zone_index]+6*i+2]/pcell->area;
1724 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T);
1725 f[zofs[zone_index]+6*i+4] = f[zofs[zone_index]+6*i+4]/pcell->area;
1726 f[zofs[zone_index]+6*i+5] = f[zofs[zone_index]+6*i+5]/pcell->area;
1727
1728 if(ODE_F.TimeDependent==true)
1729 {
1730 //second order
1731 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1732 {
1733 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1734 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1735 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1736 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1737 PetscScalar Wnt = (2-r)/(1-r)*1.5*ni*kb*Tn-1.0/(r*(1-r))*1.5*fs[i].n*kb*fs[i].Tn+(1-r)/r*1.5*fs[i].n_last*kb*fs[i].Tn_last;
1738 PetscScalar Wpt = (2-r)/(1-r)*1.5*pi*kb*Tp-1.0/(r*(1-r))*1.5*fs[i].p*kb*fs[i].Tp+(1-r)/r*1.5*fs[i].p_last*kb*fs[i].Tp_last;
1739 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1740 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1741 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1742 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1743 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
1744 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
1745 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1746 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1747 }
1748 else //first order
1749 {
1750 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1751 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1752 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1753 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1754 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
1755 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
1756 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1757 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1758 }
1759 }
1760 }
1761
1762
1763
1764 void SMCZone::F3E_ddm_homojunction(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1765 SMCZone *pz, int n)
1766 {
1767 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1768 PetscScalar kb = mt->kb;
1769 PetscScalar e = mt->e;
1770 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1771 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1772 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1773 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1774 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1775 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1776 if(zone_index > pz->pzone->zone_index)
1777 {
1778 f[zofs[zone_index]+6*i+0] = Vi - x[zofs[pz->zone_index]+6*n+0];
1779 f[zofs[zone_index]+6*i+1] = ni - x[zofs[pz->zone_index]+6*n+1];
1780 f[zofs[zone_index]+6*i+2] = pi - x[zofs[pz->zone_index]+6*n+2];
1781 f[zofs[zone_index]+6*i+3] = Ti - x[zofs[pz->zone_index]+6*n+3];
1782 f[zofs[zone_index]+6*i+4] = x[zofs[zone_index]+6*i+4] - x[zofs[pz->zone_index]+6*n+4];
1783 f[zofs[zone_index]+6*i+5] = x[zofs[zone_index]+6*i+5] - x[zofs[pz->zone_index]+6*n+5];
1784 return;
1785 }
1786
1787 //process half cell in local zone
1788 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
1789
1790 PetscScalar Na = aux[i].Na;
1791 PetscScalar Nd = aux[i].Nd;
1792 PetscScalar grad_P=0,grad_T=0;
1793 PetscScalar Fphi=0,Fn=0,Fp=0,Sn=0,Sp=0;
1794
1795 Fphi = f[zofs[zone_index]+6*i+0];
1796 Fn = f[zofs[zone_index]+6*i+1];
1797 Fp = f[zofs[zone_index]+6*i+2];
1798 grad_T=f[zofs[zone_index]+6*i+3];
1799 Sn = f[zofs[zone_index]+6*i+4];
1800 Sp = f[zofs[zone_index]+6*i+5];
1801 //process another half cell in dornor zone
1802 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1803
1804 Fphi += f[zofs[pz->zone_index]+6*n+0];
1805 Fn += f[zofs[pz->zone_index]+6*n+1];
1806 Fp += f[zofs[pz->zone_index]+6*n+2];
1807 grad_T+=f[zofs[pz->zone_index]+6*n+3];
1808 Sn += f[zofs[pz->zone_index]+6*n+4];
1809 Sp += f[zofs[pz->zone_index]+6*n+5];
1810
1811 PetscScalar area = pcell->area+ncell->area;
1812 f[zofs[zone_index]+6*i+0] = Fphi/area + e/aux[i].eps*(pi-ni+Nd-Na);
1813 f[zofs[zone_index]+6*i+1] = Fn/area;
1814 f[zofs[zone_index]+6*i+2] = Fp/area;
1815 f[zofs[zone_index]+6*i+3] = grad_T/area;
1816 f[zofs[zone_index]+6*i+4] = Sn/area;
1817 f[zofs[zone_index]+6*i+5] = Sp/area;
1818
1819 if(ODE_F.TimeDependent==true)
1820 {
1821 //second order
1822 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1823 {
1824 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1825 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1826 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1827 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1828 PetscScalar Wnt = (2-r)/(1-r)*1.5*ni*kb*Tn-1.0/(r*(1-r))*1.5*fs[i].n*kb*fs[i].Tn+(1-r)/r*1.5*fs[i].n_last*kb*fs[i].Tn_last;
1829 PetscScalar Wpt = (2-r)/(1-r)*1.5*pi*kb*Tp-1.0/(r*(1-r))*1.5*fs[i].p*kb*fs[i].Tp+(1-r)/r*1.5*fs[i].p_last*kb*fs[i].Tp_last;
1830 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1831 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1832 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1833 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1834 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1835 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1836 }
1837 else //first order
1838 {
1839 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1840 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1841 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1842 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1843 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1844 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1845 }
1846 }
1847 }
1848
1849
1850
1851
1852 void SMCZone::F3E_om_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1853 {
1854 int equ_num = 6;
1855 int size = pzone->davcell.size();
1856 int bc_index = electrode[i];
1857 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
1858 PetscScalar e = mt->e;
1859 PetscScalar current=0;
1860
1861 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1862 {
1863 int node=bc[bc_index].psegment->node_array[j];
1864 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1865 PetscScalar Vi = x[zofs[zone_index]+6*node+0]; //potential of node i
1866 current += DeviceDepth*(f[zofs[zone_index]+6*node+1]-f[zofs[zone_index]+6*node+2]);
1867 for(int k=0;k<pcell->nb_num;k++)
1868 {
1869 int nb = pcell->nb_array[k];
1870 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
1871 //displacement current
1872 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1873 }
1874 }
1875 //if the electrode connect to another electrode
1876 if(pbc->inner_connect!=-1)
1877 {
1878 int connect_to = pbc->inner_connect;
1879 int connect_zone = pbc->connect_zone;
1880 int connect_elec=-1;
1881 SMCZone * pz = dynamic_cast <SMCZone * > (pbc->pzonedata);
1882 for(int j=0;j<pz->electrode.size();j++)
1883 {
1884 if(bc[pz->electrode[j]].BCType==OhmicContact)
1885 {
1886 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pz->electrode[j]));
1887 if(om_bc->inner_connect==bc_index) connect_elec=j;
1888 }
1889 }
1890 int connect_index=zofs[connect_zone]+equ_num*pz->pzone->davcell.size()+connect_elec;
1891 if(bc_index < connect_to)
1892 {
1893 f[zofs[zone_index]+equ_num*size+i]+=x[zofs[zone_index]+equ_num*size+i]; //voltage equation
1894 f[connect_index]+=current; //current equation
1895 }
1896 else
1897 {
1898 f[zofs[zone_index]+equ_num*size+i]+=current; //current equation
1899 f[connect_index]+=-x[zofs[zone_index]+equ_num*size+i];//voltage equation
1900 }
1901 pbc->Set_Current_new(current);
1902 }
1903 //the electrode has an external voltage circult
1904 else if(pbc->electrode_type==VoltageBC)
1905 {
1906 PetscScalar V = pbc->Vapp;
1907 PetscScalar R = pbc->R;
1908 PetscScalar C = pbc->C;
1909 PetscScalar L = pbc->L;
1910 PetscScalar In = pbc->current;
1911 PetscScalar Icn = pbc->cap_current;
1912 PetscScalar Pn = pbc->potential;
1913 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1914 f[zofs[zone_index]+equ_num*size+i]=(L/ODE_F.dt+R)*current-V+(1+(L/ODE_F.dt+R)*C/ODE_F.dt)*P
1915 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1916 pbc->Set_Current_new(current);
1917 }
1918 //the electrode has an external current circult
1919 else if(pbc->electrode_type==CurrentBC)
1920 {
1921 PetscScalar I = pbc->Iapp;
1922 PetscScalar C = pbc->C;
1923 f[zofs[zone_index]+equ_num*size+i]=current - I;
1924 pbc->Set_Current_new(I);
1925 }
1926 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1927 }
1928
1929 void SMCZone::F3E_stk_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1930 {
1931 int equ_num = 6;
1932 int size = pzone->davcell.size();
1933 int bc_index = electrode[i];
1934 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
1935
1936 PetscScalar current=0;
1937 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1938 {
1939 int node = bc[bc_index].psegment->node_array[j];
1940 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1941 PetscScalar Vi = x[zofs[zone_index]+6*node+0]; //electron density of node i
1942 PetscScalar ni = x[zofs[zone_index]+6*node+1]; //electron density of node i
1943 PetscScalar pi = x[zofs[zone_index]+6*node+2]; //hole density of node i
1944 PetscScalar Ti = x[zofs[zone_index]+6*node+3];
1945 mt->mapping(&pzone->danode[node],&aux[node],ODE_F.clock);
1946 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
1947 PetscScalar Jsn = mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[node].affinity-deltaVB);
1948 PetscScalar Jsp = mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[node].affinity+deltaVB);
1949 current += -Jsn*0.5*pcell->ilen[0]*DeviceDepth;
1950 current += -Jsp*0.5*pcell->ilen[0]*DeviceDepth;
1951 current += -Jsn*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
1952 current += -Jsp*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
1953 for(int k=0;k<pcell->nb_num;k++)
1954 {
1955 int nb = pcell->nb_array[k];
1956 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
1957 //displacement current
1958 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1959 }
1960 }
1961 if(pbc->electrode_type==VoltageBC)
1962 {
1963 PetscScalar R = pbc->R;
1964 PetscScalar C = pbc->C;
1965 PetscScalar L = pbc->L;
1966 PetscScalar V = pbc->Vapp;
1967 PetscScalar In = pbc->current;
1968 PetscScalar Icn = pbc->cap_current;
1969 PetscScalar Pn = pbc->potential;
1970 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1971 f[zofs[zone_index]+equ_num*size+i]=(L/ODE_F.dt+R)*current-V+(1+(L/ODE_F.dt+R)*C/ODE_F.dt)*P
1972 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1973 pbc->Set_Current_new(current);
1974 }
1975 else if(pbc->electrode_type==CurrentBC)
1976 {
1977 PetscScalar I = pbc->Iapp;
1978 PetscScalar C = pbc->C;
1979 f[zofs[zone_index]+equ_num*size+i]=current - I;
1980 pbc->Set_Current_new(I);
1981 }
1982 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1983 }
1984
1985
1986 void SMCZone::F3E_ins_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1987 {
1988 int equ_num = 6;
1989 int size = pzone->davcell.size();
1990 int bc_index = electrode[i];
1991 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
1992
1993 PetscScalar current=0;
1994 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1995 {
1996 int node = bc[bc_index].psegment->node_array[j];
1997 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1998 PetscScalar Vi = x[zofs[zone_index]+6*node+0]; //electron density of node i
1999 for(int k=0;k<pcell->nb_num;k++)
2000 {
2001 int nb = pcell->nb_array[k];
2002 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
2003 //displacement current
2004 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
2005 }
2006 }
2007 if(pbc->electrode_type==VoltageBC)
2008 {
2009 PetscScalar R = pbc->R;
2010 PetscScalar C = pbc->C;
2011 PetscScalar L = pbc->L;
2012 PetscScalar V = pbc->Vapp;
2013 PetscScalar In = pbc->current;
2014 PetscScalar Icn = pbc->cap_current;
2015 PetscScalar Pn = pbc->potential;
2016 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
2017 f[zofs[zone_index]+equ_num*size+i]=(L/ODE_F.dt+R)*current-V+(1+(L/ODE_F.dt+R)*C/ODE_F.dt)*P
2018 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
2019 pbc->Set_Current_new(current);
2020 }
2021 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
2022 }
2023
2024
2025
2026 void SMCZone::J3E_ddm_inner(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs)
2027 {
2028 Mat6 A;
2029 PetscInt index[6],col[6];
2030 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2031 int z = zone_index;
2032 PetscScalar kb = mt->kb;
2033 PetscScalar e = mt->e;
2034 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2035 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2036 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2037 PetscScalar Ti = x[zofs[z]+6*i+3];
2038 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2039 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2040 PetscScalar area = pcell->area;
2041
2042 //--------------------------------
2043 index[0] = zofs[z]+6*i+0;
2044 index[1] = zofs[z]+6*i+1;
2045 index[2] = zofs[z]+6*i+2;
2046 index[3] = zofs[z]+6*i+3;
2047 index[4] = zofs[z]+6*i+4;
2048 index[5] = zofs[z]+6*i+5;
2049
2050 //--------------------------------
2051 for(int j=0;j<pcell->nb_num;j++)
2052 {
2053 int nb = pcell->nb_array[j];
2054 col[0] = zofs[z]+6*nb+0;
2055 col[1] = zofs[z]+6*nb+1;
2056 col[2] = zofs[z]+6*nb+2;
2057 col[3] = zofs[z]+6*nb+3;
2058 col[4] = zofs[z]+6*nb+4;
2059 col[5] = zofs[z]+6*nb+5;
2060 MatGetValues(*jtmp,6,index,6,col,A.m);
2061 A=A/area;
2062 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2063 }
2064
2065 MatGetValues(*jtmp,6,index,6,index,A.m);
2066 A=A/area;
2067 A.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
2068 A.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
2069
2070 if(ODE_F.TimeDependent==true)
2071 {
2072 //second order
2073 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2074 {
2075 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2076 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2077 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2078 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2079 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2080 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2081 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2082 }
2083 else //first order
2084 {
2085 A.m[7] += -1/ODE_F.dt;
2086 A.m[14] += -1/ODE_F.dt;
2087 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2088 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2089 A.m[28] += -1.5*kb/ODE_F.dt;
2090 A.m[35] += -1.5*kb/ODE_F.dt;
2091 }
2092 }
2093
2094 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2095 }
2096
2097
2098
2099 void SMCZone::J3E_ddm_neumannbc(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
2100 {
2101 Mat6 A;
2102 PetscInt index[6],col[6];
2103 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2104 NeumannBC *pbc = dynamic_cast <NeumannBC * >(bc.Get_pointer(pcell->bc_index-1));
2105 int z = zone_index;
2106 PetscScalar kb = mt->kb;
2107 PetscScalar e = mt->e;
2108 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2109 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2110 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2111 PetscScalar Ti = x[zofs[z]+6*i+3];
2112 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2113 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2114 PetscScalar area = pcell->area;
2115 PetscScalar d_grad_T_dTi = 0;
2116 //--------------------------------
2117 index[0] = zofs[z]+6*i+0;
2118 index[1] = zofs[z]+6*i+1;
2119 index[2] = zofs[z]+6*i+2;
2120 index[3] = zofs[z]+6*i+3;
2121 index[4] = zofs[z]+6*i+4;
2122 index[5] = zofs[z]+6*i+5;
2123
2124 //--------------------------------
2125 for(int j=0;j<pcell->nb_num;j++)
2126 {
2127 int nb = pcell->nb_array[j];
2128 col[0] = zofs[z]+6*nb+0;
2129 col[1] = zofs[z]+6*nb+1;
2130 col[2] = zofs[z]+6*nb+2;
2131 col[3] = zofs[z]+6*nb+3;
2132 col[4] = zofs[z]+6*nb+4;
2133 col[5] = zofs[z]+6*nb+5;
2134 MatGetValues(*jtmp,6,index,6,col,A.m);
2135 A=A/area;
2136 if(j==0||j==pcell->nb_num-1)
2137 {
2138 PetscScalar h = pbc->heat_transfer;
2139 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
2140 A.m[21] += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
2141 }
2142 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2143 }
2144
2145 MatGetValues(*jtmp,6,index,6,index,A.m);
2146 A=A/area;
2147 A.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
2148 A.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
2149
2150 A.m[21]+= d_grad_T_dTi;
2151
2152 if(ODE_F.TimeDependent==true)
2153 {
2154 //second order
2155 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2156 {
2157 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2158 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2159 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2160 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2161 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2162 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2163 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2164 }
2165 else //first order
2166 {
2167 A.m[7] += -1/ODE_F.dt;
2168 A.m[14] += -1/ODE_F.dt;
2169 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2170 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2171 A.m[28] += -1.5*kb/ODE_F.dt;
2172 A.m[35] += -1.5*kb/ODE_F.dt;
2173 }
2174 }
2175
2176 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2177 }
2178
2179
2180
2181
2182 void SMCZone::J3E_ddm_ombc_segment(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
2183 {
2184 Mat6 A;
2185 PetscInt index[6];
2186 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2187 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
2188 int z = zone_index;
2189 PetscScalar kb = mt->kb;
2190 PetscScalar e = mt->e;
2191 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2192 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2193 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2194 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2195 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2196 PetscScalar area = pcell->area;
2197 PetscScalar d_grad_T_dTi = 0;
2198 index[0] = zofs[z]+6*i+0;
2199 index[1] = zofs[z]+6*i+1;
2200 index[2] = zofs[z]+6*i+2;
2201 index[3] = zofs[z]+6*i+3;
2202 index[4] = zofs[z]+6*i+4;
2203 index[5] = zofs[z]+6*i+5;
2204
2205 for(int j=0;j<pcell->nb_num;j++)
2206 {
2207 int nb = pcell->nb_array[j];
2208 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2209 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
2210 //-------------------------------------
2211 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j]/area;
2212 //-------------------------------------
2213 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]/area; //dfun(3)/dT(r)
2214 if(j==0||j==pcell->nb_num-1)
2215 {
2216 PetscScalar h = pbc->heat_transfer;
2217 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
2218 d_grad_T_dTj += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
2219 }
2220 MatSetValue(*jac,zofs[z]+6*i+3,zofs[z]+6*nb+3,d_grad_T_dTj,INSERT_VALUES);
2221 }
2222
2223 Set_Mat6_I(A);
2224 A.m[21] = d_grad_T_dTi; //dfun(3)/dT(i)
2225
2226 A.m[25] = -Ti;
2227 A.m[27] = -ni;
2228
2229 A.m[32] = -Ti;
2230 A.m[33] = -pi;
2231 if(ODE_F.TimeDependent==true)
2232 {
2233 //second order
2234 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2235 {
2236 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2237 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2238 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2239 }
2240 else //first order
2241 {
2242 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2243 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2244 }
2245 }
2246 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2247
2248 int om_equ;
2249 int equ_num = 6;
2250 int size = pzone->davcell.size();
2251
2252 for(int j=0;j<electrode.size();j++)
2253 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
2254
2255 MatSetValue(*jac,zofs[zone_index]+6*i+0,zofs[zone_index]+equ_num*size+om_equ,-1,INSERT_VALUES);
2256 }
2257
2258
2259 void SMCZone::J3E_ddm_ombc_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc,
2260 ElZone *pz, int n)
2261 {
2262 Mat6 A;
2263 PetscInt index[6];
2264 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2265 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
2266 int z = zone_index;
2267 PetscScalar kb = mt->kb;
2268 PetscScalar e = mt->e;
2269 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2270 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2271 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2272 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2273 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2274 PetscScalar area = pcell->area;
2275 PetscScalar d_grad_T_dTi = 0;
2276 index[0] = zofs[z]+6*i+0;
2277 index[1] = zofs[z]+6*i+1;
2278 index[2] = zofs[z]+6*i+2;
2279 index[3] = zofs[z]+6*i+3;
2280 index[4] = zofs[z]+6*i+4;
2281 index[5] = zofs[z]+6*i+5;
2282
2283 for(int j=0;j<pcell->nb_num;j++)
2284 {
2285 int nb = pcell->nb_array[j];
2286 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2287 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
2288 //-------------------------------------
2289 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j];
2290 //-------------------------------------
2291 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]; //dfun(3)/dT(r)
2292 MatSetValue(*jac,zofs[z]+6*i+3,zofs[z]+6*nb+3,d_grad_T_dTj,INSERT_VALUES);
2293 }
2294 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2295 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
2296 for(int j=0;j<ncell->nb_num;j++)
2297 {
2298 int nb = ncell->nb_array[j];
2299 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
2300 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
2301 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
2302 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
2303 MatSetValue(*jac,zofs[z]+6*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
2304 }
2305
2306 Set_Mat6_I(A);
2307 A.m[21] = d_grad_T_dTi; //dfun(3)/dT(i)
2308
2309 A.m[25] = -Ti;
2310 A.m[27] = -ni;
2311
2312 A.m[32] = -Ti;
2313 A.m[33] = -pi;
2314
2315 if(ODE_F.TimeDependent==true)
2316 {
2317 //second order
2318 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2319 {
2320 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2321 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2322 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2323 A.m[21] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
2324 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
2325 }
2326 else //first order
2327 {
2328 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2329 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2330 A.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
2331 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
2332 }
2333 }
2334
2335 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2336
2337 int om_equ;
2338 int equ_num = 6;
2339 int size = pzone->davcell.size();
2340
2341 for(int j=0;j<electrode.size();j++)
2342 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
2343
2344 MatSetValue(*jac,zofs[zone_index]+6*i+0,zofs[zone_index]+equ_num*size+om_equ,-1,INSERT_VALUES);
2345 }
2346
2347
2348 void SMCZone::J3E_ddm_stkbc_segment(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
2349 {
2350 Mat6 A;
2351 PetscInt index[6],col[6];
2352 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2353 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
2354 int z = zone_index;
2355 int equ_num = 6;
2356 int size = pzone->davcell.size();
2357
2358 PetscScalar kb = mt->kb;
2359 PetscScalar e = mt->e;
2360 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
2361 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
2362 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2363 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2364 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2365
2366 PetscScalar area = pcell->area;
2367 PetscScalar d_grad_T_dTi = 0;
2368 PetscScalar d_Fn_dni = 0;
2369 PetscScalar d_Fn_dTi = 0;
2370 PetscScalar d_Fp_dpi = 0;
2371 PetscScalar d_Fp_dTi = 0;
2372 //--------------------------------
2373 index[0] = zofs[z]+6*i+0;
2374 index[1] = zofs[z]+6*i+1;
2375 index[2] = zofs[z]+6*i+2;
2376 index[3] = zofs[z]+6*i+3;
2377 index[4] = zofs[z]+6*i+4;
2378 index[5] = zofs[z]+6*i+5;
2379 //--------------------------------
2380 for(int j=0;j<pcell->nb_num;j++)
2381 {
2382 int nb = pcell->nb_array[j];
2383 PetscScalar nj = x[zofs[z]+6*nb+1]; //electron density of nb node
2384 PetscScalar pj = x[zofs[z]+6*nb+2]; //hole density of nb node
2385 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2386 col[0] = zofs[z]+6*nb+0;
2387 col[1] = zofs[z]+6*nb+1;
2388 col[2] = zofs[z]+6*nb+2;
2389 col[3] = zofs[z]+6*nb+3;
2390 col[4] = zofs[z]+6*nb+4;
2391 col[5] = zofs[z]+6*nb+5;
2392 MatGetValues(*jtmp,6,index,6,col,A.m);
2393 A.m[0] = 0.0;
2394 A=A/area;
2395 if(j==0||j==pcell->nb_num-1)
2396 {
2397 d_Fn_dni += mt->band->pdSchottyJsn_pdn (ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2398 d_Fn_dTi += mt->band->pdSchottyJsn_pdTl(ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2399 d_Fp_dpi += mt->band->pdSchottyJsp_pdp (pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2400 d_Fp_dTi += mt->band->pdSchottyJsp_pdTl(pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2401 PetscScalar h = pbc->heat_transfer;
2402 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
2403 A.m[21] += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
2404 }
2405 A.m[24] = 0.0;
2406 A.m[25] = 0.0;
2407 A.m[26] = 0.0;
2408 A.m[27] = 0.0;
2409 A.m[28] = 0.0;
2410 A.m[29] = 0.0;
2411
2412 A.m[30] = 0.0;
2413 A.m[31] = 0.0;
2414 A.m[32] = 0.0;
2415 A.m[33] = 0.0;
2416 A.m[34] = 0.0;
2417 A.m[35] = 0.0;
2418 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2419 }
2420
2421 MatGetValues(*jtmp,6,index,6,index,A.m);
2422 A=A/area;
2423
2424 A.m[0] = 1.0;
2425
2426 A.m[7] += d_Fn_dni; //dfun(1)/dn(i)
2427 A.m[9] += d_Fn_dTi;
2428
2429 A.m[14]+= -d_Fp_dpi; //dfun(2)/dp(i)
2430 A.m[15]+= -d_Fp_dTi;
2431
2432 A.m[21]+= d_grad_T_dTi; //dfun(3)/dT(i)
2433
2434 A.m[24] = 0;
2435 A.m[25] = -Ti;
2436 A.m[26] = 0;
2437 A.m[27] = -ni;
2438 A.m[28] = 1.0;
2439 A.m[29] = 0;
2440
2441 A.m[30] = 0.0;
2442 A.m[31] = 0.0;
2443 A.m[32] = -Ti;
2444 A.m[33] = -pi;
2445 A.m[34] = 0.0;
2446 A.m[35] = 1.0;
2447 if(ODE_F.TimeDependent==true)
2448 {
2449 //second order
2450 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2451 {
2452 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2453 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2454 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2455 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2456 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2457 }
2458 else //first order
2459 {
2460 A.m[7] += -1/ODE_F.dt;
2461 A.m[14] += -1/ODE_F.dt;
2462 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2463 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2464 }
2465 }
2466 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2467
2468 int stk_equ;
2469 for(int j=0;j<electrode.size();j++)
2470 if(electrode[j]==pcell->bc_index-1)
2471 {
2472 stk_equ=j;
2473 break;
2474 }
2475 MatSetValue(*jac,zofs[zone_index]+6*i,zofs[zone_index]+equ_num*size+stk_equ,-1.0,INSERT_VALUES);
2476 }
2477
2478
2479 void SMCZone::J3E_ddm_stkbc_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc,
2480 ElZone *pz, int n)
2481 {
2482 Mat6 A;
2483 PetscInt index[6],col[6];
2484 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2485 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
2486 int z = zone_index;
2487 int equ_num = 6;
2488 int size = pzone->davcell.size();
2489 PetscScalar area = pcell->area;
2490 Vec6 scale;
2491 Set_Vec6(scale,area,area,area,1.0,area,area);
2492
2493 PetscScalar kb = mt->kb;
2494 PetscScalar e = mt->e;
2495 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
2496 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
2497 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2498 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2499 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2500
2501 PetscScalar d_grad_T_dTi = 0;
2502 PetscScalar d_Fn_dni = 0;
2503 PetscScalar d_Fn_dTi = 0;
2504 PetscScalar d_Fp_dpi = 0;
2505 PetscScalar d_Fp_dTi = 0;
2506 //--------------------------------
2507 index[0] = zofs[z]+6*i+0;
2508 index[1] = zofs[z]+6*i+1;
2509 index[2] = zofs[z]+6*i+2;
2510 index[3] = zofs[z]+6*i+3;
2511 index[4] = zofs[z]+6*i+4;
2512 index[5] = zofs[z]+6*i+5;
2513 //--------------------------------
2514 for(int j=0;j<pcell->nb_num;j++)
2515 {
2516 int nb = pcell->nb_array[j];
2517 PetscScalar nj = x[zofs[z]+6*nb+1]; //electron density of nb node
2518 PetscScalar pj = x[zofs[z]+6*nb+2]; //hole density of nb node
2519 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2520 col[0] = zofs[z]+6*nb+0;
2521 col[1] = zofs[z]+6*nb+1;
2522 col[2] = zofs[z]+6*nb+2;
2523 col[3] = zofs[z]+6*nb+3;
2524 col[4] = zofs[z]+6*nb+4;
2525 col[5] = zofs[z]+6*nb+5;
2526 MatGetValues(*jtmp,6,index,6,col,A.m);
2527 A.m[0] = 0.0;
2528 A=A/scale;
2529 if(j==0||j==pcell->nb_num-1)
2530 {
2531 d_Fn_dni += mt->band->pdSchottyJsn_pdn (ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2532 d_Fn_dTi += mt->band->pdSchottyJsn_pdTl(ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2533 d_Fp_dpi += mt->band->pdSchottyJsp_pdp (pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2534 d_Fp_dTi += mt->band->pdSchottyJsp_pdTl(pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2535 }
2536 A.m[24] = 0.0;
2537 A.m[25] = 0.0;
2538 A.m[26] = 0.0;
2539 A.m[27] = 0.0;
2540 A.m[28] = 0.0;
2541 A.m[29] = 0.0;
2542
2543 A.m[30] = 0.0;
2544 A.m[31] = 0.0;
2545 A.m[32] = 0.0;
2546 A.m[33] = 0.0;
2547 A.m[34] = 0.0;
2548 A.m[35] = 0.0;
2549 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2550 }
2551
2552 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2553 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
2554 for(int j=0;j<ncell->nb_num;j++)
2555 {
2556 int nb = ncell->nb_array[j];
2557 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
2558 PetscScalar kapa = pz->mt->thermal->HeatConduction(fs[i].T);//(0.5*(Ti+Tj_n));
2559 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
2560 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
2561 MatSetValue(*jac,zofs[z]+6*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
2562 }
2563
2564 MatGetValues(*jtmp,6,index,6,index,A.m);
2565 A=A/scale;
2566
2567 A.m[0] = 1.0;
2568
2569 A.m[7] += d_Fn_dni; //dfun(1)/dn(i)
2570 A.m[9] += d_Fn_dTi;
2571
2572 A.m[14]+= -d_Fp_dpi; //dfun(2)/dp(i)
2573 A.m[15]+= -d_Fp_dTi;
2574
2575 A.m[21]+= d_grad_T_dTi; //dfun(3)/dT(i)
2576
2577 A.m[24] = 0;
2578 A.m[25] = -Ti;
2579 A.m[26] = 0;
2580 A.m[27] = -ni;
2581 A.m[28] = 1.0;
2582 A.m[29] = 0;
2583
2584 A.m[30] = 0.0;
2585 A.m[31] = 0.0;
2586 A.m[32] = -Ti;
2587 A.m[33] = -pi;
2588 A.m[34] = 0.0;
2589 A.m[35] = 1.0;
2590
2591 if(ODE_F.TimeDependent==true)
2592 {
2593 //second order
2594 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2595 {
2596 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2597 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2598 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2599 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2600 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2601 A.m[21] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
2602 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
2603 }
2604 else //first order
2605 {
2606 A.m[7] += -1/ODE_F.dt;
2607 A.m[14] += -1/ODE_F.dt;
2608 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2609 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2610 A.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
2611 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
2612 }
2613 }
2614
2615 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2616
2617 int stk_equ;
2618 for(int j=0;j<electrode.size();j++)
2619 if(electrode[j]==pcell->bc_index-1)
2620 {
2621 stk_equ=j;
2622 break;
2623 }
2624 MatSetValue(*jac,zofs[zone_index]+6*i,zofs[zone_index]+equ_num*size+stk_equ,-1.0,INSERT_VALUES);
2625 }
2626
2627
2628 void SMCZone::J3E_ddm_insulator_gate(int i,PetscScalar *x, Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
2629 {
2630 Mat6 A;
2631 PetscInt index[6],col[6];
2632 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2633 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
2634 int z = zone_index;
2635 int equ_num = 6;
2636 int size = pzone->davcell.size();
2637 PetscScalar kb = mt->kb;
2638 PetscScalar e = mt->e;
2639 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2640 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2641 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2642 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2643 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2644 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2645
2646 PetscScalar area = pcell->area;
2647 PetscScalar d_grad_P_dVi = 0;
2648 PetscScalar d_grad_P_dVapp = 0;
2649 PetscScalar d_grad_T_dTi = 0;
2650 int ins_equ;
2651 for(int j=0;j<electrode.size();j++)
2652 if(electrode[j]==pcell->bc_index-1)
2653 {ins_equ=j;break;}
2654 //--------------------------------
2655 index[0] = zofs[z]+6*i+0;
2656 index[1] = zofs[z]+6*i+1;
2657 index[2] = zofs[z]+6*i+2;
2658 index[3] = zofs[z]+6*i+3;
2659 index[4] = zofs[z]+6*i+4;
2660 index[5] = zofs[z]+6*i+5;
2661 //--------------------------------
2662 for(int j=0;j<pcell->nb_num;j++)
2663 {
2664 int nb = pcell->nb_array[j];
2665 PetscScalar Vj = x[zofs[z]+6*nb+0]; //potential of nb node
2666 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2667 col[0] = zofs[z]+6*nb+0;
2668 col[1] = zofs[z]+6*nb+1;
2669 col[2] = zofs[z]+6*nb+2;
2670 col[3] = zofs[z]+6*nb+3;
2671 col[4] = zofs[z]+6*nb+4;
2672 col[5] = zofs[z]+6*nb+5;
2673
2674 MatGetValues(*jtmp,6,index,6,col,A.m);
2675 A=A/area;
2676 if(j==0||j==pcell->nb_num-1)
2677 {
2678 PetscScalar Thick = pbc->Thick;
2679 PetscScalar eps_ox = mt->eps0*pbc->eps;
2680 PetscScalar eps = aux[i].eps;
2681 PetscScalar s=eps_ox/eps/Thick;
2682 d_grad_P_dVi += -0.5*pcell->ilen[j]*0.25*s*3/area;
2683 d_grad_P_dVapp += 0.5*pcell->ilen[j]*s/area;
2684 A.m[0]+= -0.5*pcell->ilen[j]*0.25*s/area;
2685 PetscScalar h = pbc->heat_transfer;
2686 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/area;
2687 A.m[21] += -0.5*pcell->ilen[j]*0.25*h/area;
2688 }
2689 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2690
2691 }
2692 //fun(0) is the poisson's equation
2693 //fun(1) is the continuous equation of electron
2694 //fun(2) is the continuous equation of hole
2695
2696 MatGetValues(*jtmp,6,index,6,index,A.m);
2697 A=A/area;
2698 A.m[0] += d_grad_P_dVi;
2699 A.m[1] += -e/aux[i].eps; //df(phi)/dn(i)
2700 A.m[2] += e/aux[i].eps; //df(phi)/dp(i)
2701
2702 A.m[21]+= d_grad_T_dTi;
2703
2704 if(ODE_F.TimeDependent==true)
2705 {
2706 //second order
2707 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2708 {
2709 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2710 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2711 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2712 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2713 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2714 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2715 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2716 }
2717 else //first order
2718 {
2719 A.m[7] += -1/ODE_F.dt;
2720 A.m[14] += -1/ODE_F.dt;
2721 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2722 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2723 A.m[28] += -1.5*kb/ODE_F.dt;
2724 A.m[35] += -1.5*kb/ODE_F.dt;
2725 }
2726 }
2727 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2728 MatSetValue(*jac,zofs[z]+6*i,zofs[z]+equ_num*size+ins_equ,d_grad_P_dVapp,INSERT_VALUES);
2729
2730 }
2731
2732
2733 void SMCZone::J3E_ddm_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
2734 ISZone *pz, int n)
2735 {
2736 Mat6 A;
2737 PetscInt index[6],col[6];
2738 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2739 InsulatorInterfaceBC *pbc = dynamic_cast<InsulatorInterfaceBC * >(bc.Get_pointer(pcell->bc_index-1));
2740 int z = zone_index;
2741 PetscScalar kb = mt->kb;
2742 PetscScalar e = mt->e;
2743 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2744 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2745 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2746 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2747 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2748 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2749
2750 PetscScalar area = pcell->area;
2751 PetscScalar d_grad_P_dVi = 0;
2752 PetscScalar d_grad_T_dTi = 0;
2753 Vec6 scale;
2754 Set_Vec6(scale,1.0,area,area,1.0,area,area);
2755 //--------------------------------
2756 index[0] = zofs[z]+6*i+0;
2757 index[1] = zofs[z]+6*i+1;
2758 index[2] = zofs[z]+6*i+2;
2759 index[3] = zofs[z]+6*i+3;
2760 index[4] = zofs[z]+6*i+4;
2761 index[5] = zofs[z]+6*i+5;
2762 //--------------------------------
2763 for(int j=0;j<pcell->nb_num;j++)
2764 {
2765 int nb = pcell->nb_array[j];
2766
2767 PetscScalar Vj = x[zofs[z]+6*nb+0]; //potential of nb node
2768 PetscScalar Tj = x[zofs[z]+6*nb+3];
2769 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
2770 d_grad_P_dVi += -aux[i].eps*pcell->elen[j]/pcell->ilen[j];
2771 col[0] = zofs[z]+6*nb+0;
2772 col[1] = zofs[z]+6*nb+1;
2773 col[2] = zofs[z]+6*nb+2;
2774 col[3] = zofs[z]+6*nb+3;
2775 col[4] = zofs[z]+6*nb+4;
2776 col[5] = zofs[z]+6*nb+5;
2777 MatGetValues(*jtmp,6,index,6,col,A.m);
2778 A=A/scale;
2779 A.m[0] = aux[i].eps*pcell->elen[j]/pcell->ilen[j]; //dfun(0)/dP(r)
2780 A.m[1] = 0; //dfun(0)/dn(r)
2781 A.m[2] = 0; //dfun(0)/dp(r)
2782 A.m[3] = 0; //dfun(0)/dT(r)
2783 A.m[4] = 0;
2784 A.m[5] = 0;
2785 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2786 }
2787
2788 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2789 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
2790 for(int j=0;j<ncell->nb_num;j++)
2791 {
2792 int nb = ncell->nb_array[j];
2793 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
2794 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
2795 d_grad_P_dVi += -pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
2796 PetscScalar d_grad_P_dVj = pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
2797 MatSetValue(*jac,zofs[z]+6*i+0,zofs[pz->pzone->zone_index]+2*nb+0,d_grad_P_dVj,INSERT_VALUES);
2798 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
2799 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
2800 MatSetValue(*jac,zofs[z]+6*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
2801 }
2802
2803 //fun(0) is the poisson's equation
2804 //fun(1) is the continuous equation of electron
2805 //fun(2) is the continuous equation of hole
2806
2807 MatGetValues(*jtmp,6,index,6,index,A.m);
2808 A=A/scale;
2809 A.m[0] = d_grad_P_dVi; //dfun(0)/dP(i)
2810 A.m[1] = -e*pcell->area; //dfun(0)/dn(i)
2811 A.m[2] = e*pcell->area; //dfun(0)/dp(i)
2812 A.m[3] = 0.0; //dfun(0)/dT(i)
2813 A.m[4] = 0.0;
2814 A.m[5] = 0.0;
2815
2816 A.m[21]+= d_grad_T_dTi;
2817
2818
2819 if(ODE_F.TimeDependent==true)
2820 {
2821 //second order
2822 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2823 {
2824 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2825 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2826 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2827 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2828 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2829 A.m[21] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
2830 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
2831 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2832 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2833 }
2834 else //first order
2835 {
2836 A.m[7] += -1/ODE_F.dt;
2837 A.m[14] += -1/ODE_F.dt;
2838 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2839 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2840 A.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
2841 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
2842 A.m[28] += -1.5*kb/ODE_F.dt;
2843 A.m[35] += -1.5*kb/ODE_F.dt;
2844 }
2845 }
2846 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2847 }
2848
2849
2850
2851 void SMCZone::J3E_ddm_homojunction(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, SMCZone *pz, int n)
2852 {
2853 Mat6 A1,A2,AJ;
2854 PetscInt index[6],indexn[6],col[6];
2855 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2856 int z = zone_index;
2857 PetscScalar kb = mt->kb;
2858 PetscScalar e = mt->e;
2859 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2860 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2861 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2862 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2863 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2864 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2865 //--------------------------------
2866 index[0] = zofs[z]+6*i+0;
2867 index[1] = zofs[z]+6*i+1;
2868 index[2] = zofs[z]+6*i+2;
2869 index[3] = zofs[z]+6*i+3;
2870 index[4] = zofs[z]+6*i+4;
2871 index[5] = zofs[z]+6*i+5;
2872 //--------------------------------
2873 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2874 indexn[0] = zofs[pz->zone_index]+6*n+0;
2875 indexn[1] = zofs[pz->zone_index]+6*n+1;
2876 indexn[2] = zofs[pz->zone_index]+6*n+2;
2877 indexn[3] = zofs[pz->zone_index]+6*n+3;
2878 indexn[4] = zofs[pz->zone_index]+6*n+4;
2879 indexn[5] = zofs[pz->zone_index]+6*n+5;
2880 //--------------------------------
2881 PetscScalar area = pcell->area + ncell->area;
2882 //--------------------------------
2883 if(zone_index > pz->pzone->zone_index)
2884 {
2885 Set_Mat6_I(AJ);
2886 MatSetValues(*jac,6,index,6,index,AJ.m,INSERT_VALUES);
2887 AJ = -1.0*AJ;
2888 MatSetValues(*jac,6,index,6,indexn,AJ.m,INSERT_VALUES);
2889 return;
2890 }
2891 //--------------------------------
2892 for(int j=0;j<pcell->nb_num;j++)
2893 {
2894 int nb = pcell->nb_array[j];
2895 col[0] = zofs[z]+6*nb+0;
2896 col[1] = zofs[z]+6*nb+1;
2897 col[2] = zofs[z]+6*nb+2;
2898 col[3] = zofs[z]+6*nb+3;
2899 col[4] = zofs[z]+6*nb+4;
2900 col[5] = zofs[z]+6*nb+5;
2901 MatGetValues(*jtmp,6,index,6,col,A1.m);
2902 //-------------------------------------
2903 A1=A1/area;
2904 MatSetValues(*jac,6,index,6,col,A1.m,INSERT_VALUES);
2905 }
2906
2907 //process another half cell in dornor zone
2908 for(int j=0;j<ncell->nb_num;j++)
2909 {
2910 int nb = ncell->nb_array[j];
2911 col[0] = zofs[pz->zone_index]+6*nb+0;
2912 col[1] = zofs[pz->zone_index]+6*nb+1;
2913 col[2] = zofs[pz->zone_index]+6*nb+2;
2914 col[3] = zofs[pz->zone_index]+6*nb+3;
2915 col[4] = zofs[pz->zone_index]+6*nb+4;
2916 col[5] = zofs[pz->zone_index]+6*nb+5;
2917 MatGetValues(*jtmp,6,indexn,6,col,A2.m);
2918 A2=A2/area;
2919 MatSetValues(*jac,6,index,6,col,A2.m,INSERT_VALUES);
2920 }
2921
2922 //fun(0) is the poisson's equation
2923 //fun(1) is the continuous equation of electron
2924 //fun(2) is the continuous equation of hole
2925 MatGetValues(*jtmp,6,index,6,index,A1.m);
2926 MatGetValues(*jtmp,6,indexn,6,indexn,A2.m);
2927 AJ=(A1+A2)/area;
2928
2929 AJ.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
2930 AJ.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
2931
2932 if(ODE_F.TimeDependent==true)
2933 {
2934 //second order
2935 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2936 {
2937 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2938 AJ.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2939 AJ.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2940 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2941 AJ.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2942 AJ.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2943 AJ.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2944 }
2945 else //first order
2946 {
2947 AJ.m[7] += -1/ODE_F.dt;
2948 AJ.m[14] += -1/ODE_F.dt;
2949 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2950 AJ.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt;
2951 AJ.m[28] += -1.5*kb/ODE_F.dt;
2952 AJ.m[35] += -1.5*kb/ODE_F.dt;
2953 }
2954 }
2955 MatSetValues(*jac,6,index,6,index,AJ.m,INSERT_VALUES);
2956 }
2957
2958
2959 void SMCZone::J3E_om_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
2960 {
2961 int equ_num = 6;
2962 int size = pzone->davcell.size();
2963 int bc_index = electrode[i];
2964 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
2965 PetscScalar A1[6],A2[6],J[6];
2966 PetscInt index[6],col[6];
2967 int matrix_row = zofs[zone_index]+equ_num*size+i;
2968 int connect_index;
2969 if(pbc->inner_connect!=-1)
2970 {
2971 int connect_to = pbc->inner_connect;
2972 int connect_zone = pbc->connect_zone;
2973 int connect_elec=-1;
2974 SMCZone * pz = dynamic_cast <SMCZone * > (pbc->pzonedata);
2975 for(int j=0;j<pz->electrode.size();j++)
2976 {
2977 if(bc[pz->electrode[j]].BCType==OhmicContact)
2978 {
2979 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pz->electrode[j]));
2980 if(om_bc->inner_connect==bc_index) connect_elec=j;
2981 }
2982 }
2983 connect_index=zofs[connect_zone]+equ_num*pz->pzone->davcell.size()+connect_elec;
2984 }
2985 PetscScalar R = pbc->R;
2986 PetscScalar C = pbc->C;
2987 PetscScalar L = pbc->L;
2988 PetscScalar kb = mt->kb;
2989 PetscScalar e = mt->e;
2990
2991 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
2992 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
2993
2994 for(int n=0;n<bc[bc_index].psegment->node_array.size();n++)
2995 {
2996 int node=bc[bc_index].psegment->node_array[n];
2997 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
2998 PetscScalar dJdisp_dVi=0,dJdisp_dVr=0;
2999 index[0] = zofs[zone_index]+6*node+0;
3000 index[1] = zofs[zone_index]+6*node+1;
3001 index[2] = zofs[zone_index]+6*node+2;
3002 index[3] = zofs[zone_index]+6*node+3;
3003 index[4] = zofs[zone_index]+6*node+4;
3004 index[5] = zofs[zone_index]+6*node+5;
3005 for(int j=0;j<pcell->nb_num;j++)
3006 {
3007 int nb = pcell->nb_array[j];
3008 col[0] = zofs[zone_index]+6*nb+0;
3009 col[1] = zofs[zone_index]+6*nb+1;
3010 col[2] = zofs[zone_index]+6*nb+2;
3011 col[3] = zofs[zone_index]+6*nb+3;
3012 col[4] = zofs[zone_index]+6*nb+4;
3013 col[5] = zofs[zone_index]+6*nb+5;
3014
3015 MatGetValues(*jtmp,1,&index[1],6,col,A1);
3016 MatGetValues(*jtmp,1,&index[2],6,col,A2);
3017
3018 //for displacement current
3019 dJdisp_dVi += aux[node].eps/pcell->ilen[j]/ODE_F.dt*pcell->elen[j];
3020 dJdisp_dVr = -aux[node].eps/pcell->ilen[j]/ODE_F.dt*pcell->elen[j];
3021 if(pbc->inner_connect!=-1)
3022 {
3023 int connect_to = pbc->inner_connect;
3024 if(bc_index < connect_to)
3025 {
3026 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3027 J[1]=(A1[1]-A2[1])*DeviceDepth;
3028 J[2]=(A1[2]-A2[2])*DeviceDepth;
3029 J[3]=(A1[3]-A2[3])*DeviceDepth;
3030 J[4]=(A1[4]-A2[4])*DeviceDepth;
3031 J[5]=(A1[5]-A2[5])*DeviceDepth;
3032 MatSetValues(*jac,1,&connect_index,6,col,J,ADD_VALUES);
3033 }
3034 else
3035 {
3036 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3037 J[1]=(A1[1]-A2[1])*DeviceDepth;
3038 J[2]=(A1[2]-A2[2])*DeviceDepth;
3039 J[3]=(A1[3]-A2[3])*DeviceDepth;
3040 J[4]=(A1[4]-A2[4])*DeviceDepth;
3041 J[5]=(A1[5]-A2[5])*DeviceDepth;
3042 MatSetValues(*jac,1,&matrix_row,6,col,J,ADD_VALUES);
3043 }
3044 }
3045 else if(pbc->electrode_type==VoltageBC)
3046 {
3047 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth*(L/ODE_F.dt+R);
3048 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
3049 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
3050 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
3051 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
3052 J[5]=(A1[5]-A2[5])*DeviceDepth*(L/ODE_F.dt+R);
3053 MatSetValues(*jac,1,&matrix_row,6,col,J,ADD_VALUES);
3054 }
3055 else if(pbc->electrode_type==CurrentBC)
3056 {
3057 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3058 J[1]=(A1[1]-A2[1])*DeviceDepth;
3059 J[2]=(A1[2]-A2[2])*DeviceDepth;
3060 J[3]=(A1[3]-A2[3])*DeviceDepth;
3061 J[4]=(A1[4]-A2[4])*DeviceDepth;
3062 J[5]=(A1[5]-A2[5])*DeviceDepth;
3063 MatSetValues(*jac,1,&matrix_row,6,col,J,ADD_VALUES);
3064 }
3065 }
3066 MatGetValues(*jtmp,1,&index[1],6,index,A1);
3067 MatGetValues(*jtmp,1,&index[2],6,index,A2);
3068 if(pbc->inner_connect!=-1)
3069 {
3070 int connect_to = pbc->inner_connect;
3071 if(bc_index < connect_to)
3072 {
3073 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3074 J[1]=(A1[1]-A2[1])*DeviceDepth;
3075 J[2]=(A1[2]-A2[2])*DeviceDepth;
3076 J[3]=(A1[3]-A2[3])*DeviceDepth;
3077 J[4]=(A1[4]-A2[4])*DeviceDepth;
3078 J[5]=(A1[5]-A2[5])*DeviceDepth;
3079 MatSetValues(*jac,1,&connect_index,6,index,J,ADD_VALUES);
3080 }
3081 else
3082 {
3083 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3084 J[1]=(A1[1]-A2[1])*DeviceDepth;
3085 J[2]=(A1[2]-A2[2])*DeviceDepth;
3086 J[3]=(A1[3]-A2[3])*DeviceDepth;
3087 J[4]=(A1[4]-A2[4])*DeviceDepth;
3088 J[5]=(A1[5]-A2[5])*DeviceDepth;
3089 MatSetValues(*jac,1,&matrix_row,6,index,J,ADD_VALUES);
3090 }
3091 }
3092 else if(pbc->electrode_type==VoltageBC)
3093 {
3094 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth*(L/ODE_F.dt+R);
3095 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
3096 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
3097 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
3098 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
3099 J[5]=(A1[5]-A2[5])*DeviceDepth*(L/ODE_F.dt+R);
3100 MatSetValues(*jac,1,&matrix_row,6,index,J,ADD_VALUES);
3101 }
3102 else if(pbc->electrode_type==CurrentBC)
3103 {
3104 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth;
3105 J[1]=(A1[1]-A2[1])*DeviceDepth;
3106 J[2]=(A1[2]-A2[2])*DeviceDepth;
3107 J[3]=(A1[3]-A2[3])*DeviceDepth;
3108 J[4]=(A1[4]-A2[4])*DeviceDepth;
3109 J[5]=(A1[5]-A2[5])*DeviceDepth;
3110 MatSetValues(*jac,1,&matrix_row,6,index,J,ADD_VALUES);
3111 }
3112 }
3113
3114 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3115 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3116 if(pbc->inner_connect!=-1)
3117 {
3118 int connect_to = pbc->inner_connect;
3119 if(bc_index < connect_to)
3120 {
3121 MatSetValue(*jac,matrix_row,matrix_row,1.0,ADD_VALUES);
3122 }
3123 else
3124 {
3125 MatSetValue(*jac,connect_index,matrix_row,-1.0,ADD_VALUES);
3126 }
3127 }
3128 else if(pbc->electrode_type==VoltageBC)
3129 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
3130
3131 }
3132
3133
3134 void SMCZone::J3E_stk_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
3135 {
3136 int equ_num = 6;
3137 int size = pzone->davcell.size();
3138 int bc_index = electrode[i];
3139 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
3140 int matrix_row = zofs[zone_index]+equ_num*size+i;
3141 PetscScalar R = pbc->R;
3142 PetscScalar C = pbc->C;
3143 PetscScalar L = pbc->L;
3144 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
3145 {
3146 int node = bc[bc_index].psegment->node_array[j];
3147 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
3148 PetscScalar VB = pbc->WorkFunction-aux[node].affinity;
3149 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
3150 PetscScalar ni = x[zofs[zone_index]+6*node+1]; //electron density of node i
3151 PetscScalar pi = x[zofs[zone_index]+6*node+2]; //hole density of node i
3152 PetscScalar dJ_dni = -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[0]
3153 -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
3154 PetscScalar dJ_dpi = -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[0]
3155 -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
3156
3157 //for displacement current
3158 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3159 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3160 for(int k=0;k<pcell->nb_num;k++)
3161 {
3162 int nb = pcell->nb_array[k];
3163 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3164 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3165 if(pbc->electrode_type==VoltageBC)
3166 {
3167 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+0,dI_dVi*(L/ODE_F.dt+R),ADD_VALUES);
3168 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*nb+0,dI_dVr*(L/ODE_F.dt+R),ADD_VALUES);
3169 }
3170 else if(pbc->electrode_type==CurrentBC)
3171 {
3172 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+0,dI_dVi,ADD_VALUES);
3173 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*nb+0,dI_dVr,ADD_VALUES);
3174 }
3175 }
3176 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3177 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3178
3179 if(pbc->electrode_type==VoltageBC)
3180 {
3181 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+1,DeviceDepth*dJ_dni*(L/ODE_F.dt+R),INSERT_VALUES);
3182 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+2,DeviceDepth*dJ_dpi*(L/ODE_F.dt+R),INSERT_VALUES);
3183 }
3184 else if(pbc->electrode_type==CurrentBC)
3185 {
3186 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+1,DeviceDepth*dJ_dni,INSERT_VALUES);
3187 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+2,DeviceDepth*dJ_dpi,INSERT_VALUES);
3188 }
3189 }
3190 if(pbc->electrode_type==VoltageBC)
3191 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
3192 }
3193
3194
3195 void SMCZone::J3E_ins_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
3196 {
3197 int equ_num = 6;
3198 int size = pzone->davcell.size();
3199 int bc_index = electrode[i];
3200 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
3201 int matrix_row = zofs[zone_index]+equ_num*size+i;
3202 PetscScalar R = pbc->R;
3203 PetscScalar C = pbc->C;
3204 PetscScalar L = pbc->L;
3205
3206 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3207 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3208 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
3209 {
3210 int node = bc[bc_index].psegment->node_array[j];
3211 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
3212
3213 //for displacement current
3214
3215 for(int k=0;k<pcell->nb_num;k++)
3216 {
3217 int nb = pcell->nb_array[k];
3218 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3219 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3220 if(pbc->electrode_type==VoltageBC)
3221 {
3222 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+0,dI_dVi*(L/ODE_F.dt+R),ADD_VALUES);
3223 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*nb+0,dI_dVr*(L/ODE_F.dt+R),ADD_VALUES);
3224 }
3225 }
3226
3227 }
3228 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3229 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3230
3231 if(pbc->electrode_type==VoltageBC)
3232 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
3233 }
3234
3235
3236 void SMCZone::F3E_efield_update(PetscScalar *x,vector<int> & zofs, DABC &bc, vector<BZoneData *>zonedata)
3237 {
3238 PetscScalar P_x=0,P_y=0,w=0;
3239 //calculate Ex Ey with least-squares gradient construction
3240 for(int i=0; i<pzone->davcell.size();i++)
3241 {
3242 P_x=0,P_y=0,w=0;
3243 VoronoiCell *pcell = pzone->davcell.GetPointer(i);
3244 for(int k=0;k<pcell->nb_num;k++)
3245 {
3246 const VoronoiCell *ncell = pzone->davcell.GetPointer(pcell->nb_array[k]);
3247 PetscScalar dx = ncell->x - pcell->x;
3248 PetscScalar dy = ncell->y - pcell->y;
3249 PetscScalar dP = x[zofs[zone_index]+6*pcell->nb_array[k]+0] - x[zofs[zone_index]+6*i+0];
3250 w=1.0/sqrt(dx*dx+dy*dy);
3251 P_x+=w*w*dP*dx;
3252 P_y+=w*w*dP*dy;
3253 }
3254 if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
3255 {
3256 InsulatorInterfaceBC *pbc;
3257 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
3258 int n_zone = pbc->pinterface->Find_neighbor_zone_index(zone_index);
3259 int n_node = pbc->pinterface->Find_neighbor_node_index(zone_index,i);
3260 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
3261 const VoronoiCell* dcell = pz->pzone->davcell.GetPointer(n_node);
3262 for(int k=1;k<dcell->nb_num-1;k++)
3263 {
3264 const VoronoiCell *ncell = pz->pzone->davcell.GetPointer(dcell->nb_array[k]);
3265 PetscScalar dx = ncell->x - dcell->x;
3266 PetscScalar dy = ncell->y - dcell->y;
3267 PetscScalar dP = x[zofs[n_zone]+2*dcell->nb_array[k]]- x[zofs[n_zone]+2*n_node];
3268 w=1.0/sqrt(dx*dx+dy*dy);
3269 P_x+=w*w*dP*dx;
3270 P_y+=w*w*dP*dy;
3271 }
3272 }
3273 aux[i].Ex = -(pcell->sc*P_x-pcell->sb*P_y);
3274 aux[i].Ey = -(pcell->sa*P_y-pcell->sb*P_x);
3275 }
3276 }
General-purpose Semiconductor Simulator