?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
initial checkin, based on GSS 0.46 CVS
[gss-tcad.git] / src / solver / ebm3e / semiequ3e.cc
blob63a70509f553d5a14248fb98035b81795ebbb287
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 electron_density,hole_density;
1295 int om_equ;
1296 for(int j=0;j<electrode.size();j++)
1297 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
1298
1299 f[zofs[z]+6*i+0] = Vi - kb*fs[i].T/mt->e*asinh((Nd-Na)/(2*nie)) + mt->kb*fs[i].T/mt->e*log(Nc/nie)
1300 + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T)
1301 - x[zofs[z]+equ_num*size+om_equ];
1302
1303 if(Na>Nd) //p-type
1304 {
1305 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1306 electron_density = nie*nie/hole_density;
1307 }
1308 else //n-type
1309 {
1310 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1311 hole_density = nie*nie/electron_density;
1312 }
1313 f[zofs[z]+6*i+1] = x[zofs[z]+6*i+1] - electron_density; //electron density
1314 f[zofs[z]+6*i+2] = x[zofs[z]+6*i+2] - hole_density; //hole density
1315
1316 PetscScalar grad_T=0;
1317 for(int j=0;j<pcell->nb_num;j++)
1318 {
1319 int nb = pcell->nb_array[j];
1320 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1321 PetscScalar T_mid = 0.5*(Tj+Ti);
1322 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
1323 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
1324 grad_T += kapa*pcell->elen[j]*dTdx_mid/pcell->area;
1325 if(j==0||j==pcell->nb_num-1)
1326 {
1327 PetscScalar h = pbc->heat_transfer;
1328 PetscScalar r = h*pbc->T_external;
1329 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj))/pcell->area;
1330 }
1331 }
1332 f[zofs[zone_index]+6*i+3] = grad_T;
1333 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1334 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1335
1336 if(ODE_F.TimeDependent==true)
1337 {
1338 //second order
1339 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1340 {
1341 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1342 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1343 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1344 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1345 }
1346 else //first order
1347 {
1348 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1349 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1350 }
1351 }
1352 }
1353
1354
1355 void SMCZone::F3E_ddm_ombc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1356 ElZone *pz, int n)
1357 {
1358 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1359 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
1360 int z = zone_index;
1361 int equ_num = 6;
1362 int size = pzone->davcell.size();
1363 PetscScalar kb = mt->kb;
1364 PetscScalar Na = aux[i].Na;
1365 PetscScalar Nd = aux[i].Nd;
1366 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1367 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1368 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1369 PetscScalar Ti = x[zofs[z]+6*i+3];
1370 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1371 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1372 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
1373 PetscScalar nie = mt->band->nie(fs[i].T);
1374 PetscScalar Nc = mt->band->Nc(fs[i].T);
1375 PetscScalar electron_density,hole_density;
1376 int om_equ;
1377 for(int j=0;j<electrode.size();j++)
1378 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
1379
1380 f[zofs[z]+6*i+0] = Vi - kb*fs[i].T/mt->e*asinh((Nd-Na)/(2*nie)) + mt->kb*fs[i].T/mt->e*log(Nc/nie)
1381 + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T)
1382 - x[zofs[z]+equ_num*size+om_equ];
1383
1384 if(Na>Nd) //p-type
1385 {
1386 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1387 electron_density = nie*nie/hole_density;
1388 }
1389 else //n-type
1390 {
1391 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
1392 hole_density = nie*nie/electron_density;
1393 }
1394 f[zofs[z]+6*i+1] = x[zofs[z]+6*i+1] - electron_density; //electron density
1395 f[zofs[z]+6*i+2] = x[zofs[z]+6*i+2] - hole_density; //hole density
1396
1397 PetscScalar grad_T=0;
1398 for(int j=0;j<pcell->nb_num;j++)
1399 {
1400 int nb = pcell->nb_array[j];
1401 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1402 PetscScalar T_mid = 0.5*(Tj+Ti);
1403 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
1404 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
1405 grad_T += kapa*pcell->elen[j]*dTdx_mid;
1406 }
1407 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1408 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
1409 for(int j=0;j<ncell->nb_num;j++)
1410 {
1411 int nb = ncell->nb_array[j];
1412 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
1413 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
1414 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
1415 }
1416
1417 f[zofs[zone_index]+6*i+3] = grad_T;
1418 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1419 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1420
1421
1422 if(ODE_F.TimeDependent==true)
1423 {
1424 //second order
1425 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1426 {
1427 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1428 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1429 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1430 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1431 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
1432 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
1433 }
1434 else //first order
1435 {
1436 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1437 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1438 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
1439 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
1440 }
1441 }
1442
1443 }
1444
1445
1446
1447 void SMCZone::F3E_ddm_stkbc_segment(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc)
1448 {
1449 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1450 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
1451 int z = zone_index;
1452 int equ_num = 6;
1453 int size = pzone->davcell.size();
1454 PetscScalar kb = mt->kb;
1455 PetscScalar e = mt->e;
1456 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1457 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1458 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1459 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
1460 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1461 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1462
1463 int stk_equ;
1464 for(int j=0;j<electrode.size();j++)
1465 if(electrode[j]==pcell->bc_index-1)
1466 {
1467 stk_equ=j;
1468 break;
1469 }
1470 //Schotty Barrier Lowerring
1471 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
1472 //potential
1473 f[zofs[z]+6*i+0] = x[zofs[z]+6*i+0] + pbc->WorkFunction - deltaVB - x[zofs[z]+equ_num*size+stk_equ];
1474
1475 PetscScalar Fn=0,Fp=0,grad_T=0;
1476 for(int j=0;j<pcell->nb_num;j++)
1477 {
1478 int nb = pcell->nb_array[j];
1479 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1480 if(j==0||j==pcell->nb_num-1)
1481 {
1482 Fn += mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[i].affinity - deltaVB)*0.5*pcell->ilen[j];
1483 Fp += mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[i].affinity + deltaVB)*0.5*pcell->ilen[j];
1484 PetscScalar h = pbc->heat_transfer;
1485 PetscScalar r = h*pbc->T_external;
1486 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj));
1487 }
1488 }
1489
1490 f[zofs[zone_index]+6*i+1] = (f[zofs[zone_index]+6*i+1]+Fn)/pcell->area;
1491 f[zofs[zone_index]+6*i+2] = (f[zofs[zone_index]+6*i+2]-Fp)/pcell->area;
1492 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T)/pcell->area;
1493 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1494 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1495
1496 if(ODE_F.TimeDependent==true)
1497 {
1498 //second order
1499 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1500 {
1501 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1502 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1503 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1504 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1505 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1506 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1507 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1508 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1509 }
1510 else //first order
1511 {
1512 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1513 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1514 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1515 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1516 }
1517 }
1518 }
1519
1520
1521 void SMCZone::F3E_ddm_stkbc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc,
1522 ElZone *pz, int n)
1523 {
1524 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1525 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
1526 int z = zone_index;
1527 int equ_num = 6;
1528 int size = pzone->davcell.size();
1529 PetscScalar kb = mt->kb;
1530 PetscScalar e = mt->e;
1531 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
1532 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
1533 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
1534 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
1535 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
1536 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
1537
1538 int stk_equ;
1539 for(int j=0;j<electrode.size();j++)
1540 if(electrode[j]==pcell->bc_index-1)
1541 {
1542 stk_equ=j;
1543 break;
1544 }
1545 //Schotty Barrier Lowerring
1546 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
1547 //potential
1548 f[zofs[z]+6*i+0] = x[zofs[z]+6*i+0] + pbc->WorkFunction - deltaVB - x[zofs[z]+equ_num*size+stk_equ];
1549
1550 PetscScalar Fn=0,Fp=0,grad_T=0;
1551 for(int j=0;j<pcell->nb_num;j++)
1552 {
1553 int nb = pcell->nb_array[j];
1554 if(j==0||j==pcell->nb_num-1)
1555 {
1556 Fn += mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[i].affinity - deltaVB)*0.5*pcell->ilen[j];
1557 Fp += mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[i].affinity + deltaVB)*0.5*pcell->ilen[j];
1558 }
1559 }
1560
1561 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1562 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
1563 for(int j=0;j<ncell->nb_num;j++)
1564 {
1565 int nb = ncell->nb_array[j];
1566 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
1567 PetscScalar kapa = pz->mt->thermal->HeatConduction(fs[i].T);//(0.5*(Ti+Tj_n));
1568 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
1569 }
1570
1571 f[zofs[zone_index]+6*i+1] = (f[zofs[zone_index]+6*i+1]+Fn)/pcell->area;
1572 f[zofs[zone_index]+6*i+2] = (f[zofs[zone_index]+6*i+2]-Fp)/pcell->area;
1573 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T);
1574 f[zofs[zone_index]+6*i+4] = ni*(Tn - Ti);
1575 f[zofs[zone_index]+6*i+5] = pi*(Tp - Ti);
1576
1577 if(ODE_F.TimeDependent==true)
1578 {
1579 //second order
1580 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1581 {
1582 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1583 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1584 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1585 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1586 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1587 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1588 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1589 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1590 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
1591 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
1592 }
1593 else //first order
1594 {
1595 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1596 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1597 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1598 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1599 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
1600 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
1601 }
1602 }
1603
1604 }
1605
1606
1607 void SMCZone::F3E_ddm_insulator_gate(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
1608 {
1609 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1610 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
1611 int equ_num = 6;
1612 int size = pzone->davcell.size();
1613 PetscScalar kb = mt->kb;
1614 PetscScalar e = mt->e;
1615 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1616 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1617 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1618 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1619 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1620 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1621 PetscScalar grad_P = 0, grad_T = 0;
1622 int ins_equ;
1623 for(int j=0;j<electrode.size();j++)
1624 if(electrode[j]==pcell->bc_index-1)
1625 {ins_equ=j;break;}
1626 for(int j=0;j<pcell->nb_num;j++)
1627 {
1628 int nb = pcell->nb_array[j];
1629 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
1630 PetscScalar Tj = x[zofs[zone_index]+6*nb+3]; //lattice temperature of nb node
1631
1632 if(j==0||j==pcell->nb_num-1)
1633 {
1634 PetscScalar vgate = x[zofs[zone_index]+equ_num*size+ins_equ] - pbc->WorkFunction;
1635 PetscScalar q = e*pbc->QF; //sigma is the surface change density
1636 PetscScalar Thick = pbc->Thick;
1637 PetscScalar eps_ox = mt->eps0*pbc->eps;
1638 PetscScalar eps = aux[i].eps;
1639 PetscScalar r=q/eps + eps_ox/eps/Thick*vgate;
1640 PetscScalar s=eps_ox/eps/Thick;
1641 grad_P += 0.5*pcell->ilen[j]*(r-0.25*s*(3*Vi+Vj));
1642
1643 PetscScalar h = pbc->heat_transfer;
1644 grad_T += 0.5*pcell->ilen[j]*(h*pbc->T_external-0.25*h*(3*Ti+Tj));
1645 }
1646 }
1647
1648 f[zofs[zone_index]+6*i+0] = (f[zofs[zone_index]+6*i+0]+grad_P)/pcell->area
1649 + e/aux[i].eps*((pi-aux[i].Na)+(aux[i].Nd-ni));
1650 f[zofs[zone_index]+6*i+1] = f[zofs[zone_index]+6*i+1]/pcell->area;
1651 f[zofs[zone_index]+6*i+2] = f[zofs[zone_index]+6*i+2]/pcell->area;
1652 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T)/pcell->area;
1653 f[zofs[zone_index]+6*i+4] = f[zofs[zone_index]+6*i+4]/pcell->area;
1654 f[zofs[zone_index]+6*i+5] = f[zofs[zone_index]+6*i+5]/pcell->area;
1655
1656 if(ODE_F.TimeDependent==true)
1657 {
1658 //second order
1659 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1660 {
1661 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1662 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1663 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1664 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1665 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;
1666 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;
1667 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1668 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1669 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1670 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1671 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1672 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1673 }
1674 else //first order
1675 {
1676 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1677 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1678 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1679 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1680 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1681 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1682 }
1683 }
1684 }
1685
1686
1687
1688 void SMCZone::F3E_ddm_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1689 ISZone *pz, int n)
1690 {
1691 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1692 InsulatorInterfaceBC *pbc = dynamic_cast<InsulatorInterfaceBC * >(bc.Get_pointer(pcell->bc_index-1));
1693 PetscScalar kb = mt->kb;
1694 PetscScalar e = mt->e;
1695 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1696 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1697 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1698 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1699 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1700 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1701
1702 PetscScalar grad_P = 0,grad_T=0;
1703 //semiconductor-insulator interface
1704 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1705 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
1706 for(int j=0;j<ncell->nb_num;j++)
1707 {
1708 int nb = ncell->nb_array[j];
1709 PetscScalar Vj_n = x[zofs[pz->pzone->zone_index]+2*nb+0]; //potential of nb node
1710 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
1711 grad_P += pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j]*(Vj_n-Vi);
1712 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
1713 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
1714 }
1715 PetscScalar L = (pcell->ilen[0]+pcell->ilen[pcell->nb_num-1])/2.0;
1716
1717
1718 f[zofs[zone_index]+6*i+0] = (aux[i].eps*f[zofs[zone_index]+6*i+0]+grad_P)
1719 + e*((pi-aux[i].Na)+(aux[i].Nd-ni))*pcell->area + pbc->QF*L;
1720 f[zofs[zone_index]+6*i+1] = f[zofs[zone_index]+6*i+1]/pcell->area;
1721 f[zofs[zone_index]+6*i+2] = f[zofs[zone_index]+6*i+2]/pcell->area;
1722 f[zofs[zone_index]+6*i+3] = (f[zofs[zone_index]+6*i+3]+grad_T);
1723 f[zofs[zone_index]+6*i+4] = f[zofs[zone_index]+6*i+4]/pcell->area;
1724 f[zofs[zone_index]+6*i+5] = f[zofs[zone_index]+6*i+5]/pcell->area;
1725
1726 if(ODE_F.TimeDependent==true)
1727 {
1728 //second order
1729 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1730 {
1731 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1732 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1733 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1734 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1735 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;
1736 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;
1737 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1738 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1739 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1740 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1741 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
1742 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
1743 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1744 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1745 }
1746 else //first order
1747 {
1748 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1749 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1750 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
1751 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
1752 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
1753 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
1754 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1755 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1756 }
1757 }
1758 }
1759
1760
1761
1762 void SMCZone::F3E_ddm_homojunction(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1763 SMCZone *pz, int n)
1764 {
1765 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1766 PetscScalar kb = mt->kb;
1767 PetscScalar e = mt->e;
1768 PetscScalar Vi = x[zofs[zone_index]+6*i+0]; //potential of node i
1769 PetscScalar ni = x[zofs[zone_index]+6*i+1]; //electron density of node i
1770 PetscScalar pi = x[zofs[zone_index]+6*i+2]; //hole density of node i
1771 PetscScalar Ti = x[zofs[zone_index]+6*i+3]; //lattice temperature of node i
1772 PetscScalar Tn = x[zofs[zone_index]+6*i+4]/ni;
1773 PetscScalar Tp = x[zofs[zone_index]+6*i+5]/pi;
1774 if(zone_index > pz->pzone->zone_index)
1775 {
1776 f[zofs[zone_index]+6*i+0] = Vi - x[zofs[pz->zone_index]+6*n+0];
1777 f[zofs[zone_index]+6*i+1] = ni - x[zofs[pz->zone_index]+6*n+1];
1778 f[zofs[zone_index]+6*i+2] = pi - x[zofs[pz->zone_index]+6*n+2];
1779 f[zofs[zone_index]+6*i+3] = Ti - x[zofs[pz->zone_index]+6*n+3];
1780 f[zofs[zone_index]+6*i+4] = x[zofs[zone_index]+6*i+4] - x[zofs[pz->zone_index]+6*n+4];
1781 f[zofs[zone_index]+6*i+5] = x[zofs[zone_index]+6*i+5] - x[zofs[pz->zone_index]+6*n+5];
1782 return;
1783 }
1784
1785 //process half cell in local zone
1786 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
1787
1788 PetscScalar Na = aux[i].Na;
1789 PetscScalar Nd = aux[i].Nd;
1790 PetscScalar grad_P=0,grad_T=0;
1791 PetscScalar Fphi=0,Fn=0,Fp=0,Sn=0,Sp=0;
1792
1793 Fphi = f[zofs[zone_index]+6*i+0];
1794 Fn = f[zofs[zone_index]+6*i+1];
1795 Fp = f[zofs[zone_index]+6*i+2];
1796 grad_T=f[zofs[zone_index]+6*i+3];
1797 Sn = f[zofs[zone_index]+6*i+4];
1798 Sp = f[zofs[zone_index]+6*i+5];
1799 //process another half cell in dornor zone
1800 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1801
1802 Fphi += f[zofs[pz->zone_index]+6*n+0];
1803 Fn += f[zofs[pz->zone_index]+6*n+1];
1804 Fp += f[zofs[pz->zone_index]+6*n+2];
1805 grad_T+=f[zofs[pz->zone_index]+6*n+3];
1806 Sn += f[zofs[pz->zone_index]+6*n+4];
1807 Sp += f[zofs[pz->zone_index]+6*n+5];
1808
1809 PetscScalar area = pcell->area+ncell->area;
1810 f[zofs[zone_index]+6*i+0] = Fphi/area + e/aux[i].eps*(pi-ni+Nd-Na);
1811 f[zofs[zone_index]+6*i+1] = Fn/area;
1812 f[zofs[zone_index]+6*i+2] = Fp/area;
1813 f[zofs[zone_index]+6*i+3] = grad_T/area;
1814 f[zofs[zone_index]+6*i+4] = Sn/area;
1815 f[zofs[zone_index]+6*i+5] = Sp/area;
1816
1817 if(ODE_F.TimeDependent==true)
1818 {
1819 //second order
1820 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1821 {
1822 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1823 PetscScalar nt = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1824 PetscScalar pt = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1825 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
1826 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;
1827 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;
1828 f[zofs[zone_index]+6*i+1] += -nt/(ODE_F.dt_last+ODE_F.dt);
1829 f[zofs[zone_index]+6*i+2] += -pt/(ODE_F.dt_last+ODE_F.dt);
1830 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1831 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
1832 f[zofs[zone_index]+6*i+4] += -Wnt/(ODE_F.dt_last+ODE_F.dt);
1833 f[zofs[zone_index]+6*i+5] += -Wpt/(ODE_F.dt_last+ODE_F.dt);
1834 }
1835 else //first order
1836 {
1837 f[zofs[zone_index]+6*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1838 f[zofs[zone_index]+6*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1839 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
1840 f[zofs[zone_index]+6*i+3] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
1841 f[zofs[zone_index]+6*i+4] += -(1.5*ni*kb*Tn - 1.5*fs[i].n*kb*fs[i].Tn)/ODE_F.dt;
1842 f[zofs[zone_index]+6*i+5] += -(1.5*pi*kb*Tp - 1.5*fs[i].p*kb*fs[i].Tp)/ODE_F.dt;
1843 }
1844 }
1845 }
1846
1847
1848
1849
1850 void SMCZone::F3E_om_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1851 {
1852 int equ_num = 6;
1853 int size = pzone->davcell.size();
1854 int bc_index = electrode[i];
1855 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
1856 PetscScalar e = mt->e;
1857 PetscScalar current=0;
1858
1859 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1860 {
1861 int node=bc[bc_index].psegment->node_array[j];
1862 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1863 PetscScalar Vi = x[zofs[zone_index]+6*node+0]; //potential of node i
1864 current += DeviceDepth*(f[zofs[zone_index]+6*node+1]-f[zofs[zone_index]+6*node+2]);
1865 for(int k=0;k<pcell->nb_num;k++)
1866 {
1867 int nb = pcell->nb_array[k];
1868 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
1869 //displacement current
1870 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1871 }
1872 }
1873 //if the electrode connect to another electrode
1874 if(pbc->inner_connect!=-1)
1875 {
1876 int connect_to = pbc->inner_connect;
1877 int connect_zone = pbc->connect_zone;
1878 int connect_elec=-1;
1879 SMCZone * pz = dynamic_cast <SMCZone * > (pbc->pzonedata);
1880 for(int j=0;j<pz->electrode.size();j++)
1881 {
1882 if(bc[pz->electrode[j]].BCType==OhmicContact)
1883 {
1884 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pz->electrode[j]));
1885 if(om_bc->inner_connect==bc_index) connect_elec=j;
1886 }
1887 }
1888 int connect_index=zofs[connect_zone]+equ_num*pz->pzone->davcell.size()+connect_elec;
1889 if(bc_index < connect_to)
1890 {
1891 f[zofs[zone_index]+equ_num*size+i]+=x[zofs[zone_index]+equ_num*size+i]; //voltage equation
1892 f[connect_index]+=current; //current equation
1893 }
1894 else
1895 {
1896 f[zofs[zone_index]+equ_num*size+i]+=current; //current equation
1897 f[connect_index]+=-x[zofs[zone_index]+equ_num*size+i];//voltage equation
1898 }
1899 pbc->Set_Current_new(current);
1900 }
1901 //the electrode has an external voltage circult
1902 else if(pbc->electrode_type==VoltageBC)
1903 {
1904 PetscScalar V = pbc->Vapp;
1905 PetscScalar R = pbc->R;
1906 PetscScalar C = pbc->C;
1907 PetscScalar L = pbc->L;
1908 PetscScalar In = pbc->current;
1909 PetscScalar Icn = pbc->cap_current;
1910 PetscScalar Pn = pbc->potential;
1911 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1912 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
1913 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1914 pbc->Set_Current_new(current);
1915 }
1916 //the electrode has an external current circult
1917 else if(pbc->electrode_type==CurrentBC)
1918 {
1919 PetscScalar I = pbc->Iapp;
1920 PetscScalar C = pbc->C;
1921 f[zofs[zone_index]+equ_num*size+i]=current - I;
1922 pbc->Set_Current_new(I);
1923 }
1924 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1925 }
1926
1927 void SMCZone::F3E_stk_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1928 {
1929 int equ_num = 6;
1930 int size = pzone->davcell.size();
1931 int bc_index = electrode[i];
1932 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
1933
1934 PetscScalar current=0;
1935 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1936 {
1937 int node = bc[bc_index].psegment->node_array[j];
1938 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1939 PetscScalar Vi = x[zofs[zone_index]+6*node+0]; //electron density of node i
1940 PetscScalar ni = x[zofs[zone_index]+6*node+1]; //electron density of node i
1941 PetscScalar pi = x[zofs[zone_index]+6*node+2]; //hole density of node i
1942 PetscScalar Ti = x[zofs[zone_index]+6*node+3];
1943 mt->mapping(&pzone->danode[node],&aux[node],ODE_F.clock);
1944 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
1945 PetscScalar Jsn = mt->band->SchottyJsn(ni,Ti,pbc->WorkFunction-aux[node].affinity-deltaVB);
1946 PetscScalar Jsp = mt->band->SchottyJsp(pi,Ti,pbc->WorkFunction-aux[node].affinity+deltaVB);
1947 current += -Jsn*0.5*pcell->ilen[0]*DeviceDepth;
1948 current += -Jsp*0.5*pcell->ilen[0]*DeviceDepth;
1949 current += -Jsn*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
1950 current += -Jsp*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
1951 for(int k=0;k<pcell->nb_num;k++)
1952 {
1953 int nb = pcell->nb_array[k];
1954 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
1955 //displacement current
1956 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1957 }
1958 }
1959 if(pbc->electrode_type==VoltageBC)
1960 {
1961 PetscScalar R = pbc->R;
1962 PetscScalar C = pbc->C;
1963 PetscScalar L = pbc->L;
1964 PetscScalar V = pbc->Vapp;
1965 PetscScalar In = pbc->current;
1966 PetscScalar Icn = pbc->cap_current;
1967 PetscScalar Pn = pbc->potential;
1968 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1969 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
1970 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1971 pbc->Set_Current_new(current);
1972 }
1973 else if(pbc->electrode_type==CurrentBC)
1974 {
1975 PetscScalar I = pbc->Iapp;
1976 PetscScalar C = pbc->C;
1977 f[zofs[zone_index]+equ_num*size+i]=current - I;
1978 pbc->Set_Current_new(I);
1979 }
1980 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1981 }
1982
1983
1984 void SMCZone::F3E_ins_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1985 {
1986 int equ_num = 6;
1987 int size = pzone->davcell.size();
1988 int bc_index = electrode[i];
1989 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
1990
1991 PetscScalar current=0;
1992 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1993 {
1994 int node = bc[bc_index].psegment->node_array[j];
1995 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1996 PetscScalar Vi = x[zofs[zone_index]+6*node+0]; //electron density of node i
1997 for(int k=0;k<pcell->nb_num;k++)
1998 {
1999 int nb = pcell->nb_array[k];
2000 PetscScalar Vj = x[zofs[zone_index]+6*nb+0]; //potential of nb node
2001 //displacement current
2002 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
2003 }
2004 }
2005 if(pbc->electrode_type==VoltageBC)
2006 {
2007 PetscScalar R = pbc->R;
2008 PetscScalar C = pbc->C;
2009 PetscScalar L = pbc->L;
2010 PetscScalar V = pbc->Vapp;
2011 PetscScalar In = pbc->current;
2012 PetscScalar Icn = pbc->cap_current;
2013 PetscScalar Pn = pbc->potential;
2014 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
2015 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
2016 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
2017 pbc->Set_Current_new(current);
2018 }
2019 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
2020 }
2021
2022
2023
2024 void SMCZone::J3E_ddm_inner(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs)
2025 {
2026 Mat6 A;
2027 PetscInt index[6],col[6];
2028 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2029 int z = zone_index;
2030 PetscScalar kb = mt->kb;
2031 PetscScalar e = mt->e;
2032 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2033 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2034 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2035 PetscScalar Ti = x[zofs[z]+6*i+3];
2036 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2037 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2038 PetscScalar area = pcell->area;
2039
2040 //--------------------------------
2041 index[0] = zofs[z]+6*i+0;
2042 index[1] = zofs[z]+6*i+1;
2043 index[2] = zofs[z]+6*i+2;
2044 index[3] = zofs[z]+6*i+3;
2045 index[4] = zofs[z]+6*i+4;
2046 index[5] = zofs[z]+6*i+5;
2047
2048 //--------------------------------
2049 for(int j=0;j<pcell->nb_num;j++)
2050 {
2051 int nb = pcell->nb_array[j];
2052 col[0] = zofs[z]+6*nb+0;
2053 col[1] = zofs[z]+6*nb+1;
2054 col[2] = zofs[z]+6*nb+2;
2055 col[3] = zofs[z]+6*nb+3;
2056 col[4] = zofs[z]+6*nb+4;
2057 col[5] = zofs[z]+6*nb+5;
2058 MatGetValues(*jtmp,6,index,6,col,A.m);
2059 A=A/area;
2060 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2061 }
2062
2063 MatGetValues(*jtmp,6,index,6,index,A.m);
2064 A=A/area;
2065 A.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
2066 A.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
2067
2068 if(ODE_F.TimeDependent==true)
2069 {
2070 //second order
2071 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2072 {
2073 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2074 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2075 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2076 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2077 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2078 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2079 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2080 }
2081 else //first order
2082 {
2083 A.m[7] += -1/ODE_F.dt;
2084 A.m[14] += -1/ODE_F.dt;
2085 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2086 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2087 A.m[28] += -1.5*kb/ODE_F.dt;
2088 A.m[35] += -1.5*kb/ODE_F.dt;
2089 }
2090 }
2091
2092 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2093 }
2094
2095
2096
2097 void SMCZone::J3E_ddm_neumannbc(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
2098 {
2099 Mat6 A;
2100 PetscInt index[6],col[6];
2101 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2102 NeumannBC *pbc = dynamic_cast <NeumannBC * >(bc.Get_pointer(pcell->bc_index-1));
2103 int z = zone_index;
2104 PetscScalar kb = mt->kb;
2105 PetscScalar e = mt->e;
2106 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2107 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2108 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2109 PetscScalar Ti = x[zofs[z]+6*i+3];
2110 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2111 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2112 PetscScalar area = pcell->area;
2113 PetscScalar d_grad_T_dTi = 0;
2114 //--------------------------------
2115 index[0] = zofs[z]+6*i+0;
2116 index[1] = zofs[z]+6*i+1;
2117 index[2] = zofs[z]+6*i+2;
2118 index[3] = zofs[z]+6*i+3;
2119 index[4] = zofs[z]+6*i+4;
2120 index[5] = zofs[z]+6*i+5;
2121
2122 //--------------------------------
2123 for(int j=0;j<pcell->nb_num;j++)
2124 {
2125 int nb = pcell->nb_array[j];
2126 col[0] = zofs[z]+6*nb+0;
2127 col[1] = zofs[z]+6*nb+1;
2128 col[2] = zofs[z]+6*nb+2;
2129 col[3] = zofs[z]+6*nb+3;
2130 col[4] = zofs[z]+6*nb+4;
2131 col[5] = zofs[z]+6*nb+5;
2132 MatGetValues(*jtmp,6,index,6,col,A.m);
2133 A=A/area;
2134 if(j==0||j==pcell->nb_num-1)
2135 {
2136 PetscScalar h = pbc->heat_transfer;
2137 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
2138 A.m[21] += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
2139 }
2140 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2141 }
2142
2143 MatGetValues(*jtmp,6,index,6,index,A.m);
2144 A=A/area;
2145 A.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
2146 A.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
2147
2148 A.m[21]+= d_grad_T_dTi;
2149
2150 if(ODE_F.TimeDependent==true)
2151 {
2152 //second order
2153 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2154 {
2155 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2156 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2157 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2158 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2159 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2160 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2161 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2162 }
2163 else //first order
2164 {
2165 A.m[7] += -1/ODE_F.dt;
2166 A.m[14] += -1/ODE_F.dt;
2167 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2168 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2169 A.m[28] += -1.5*kb/ODE_F.dt;
2170 A.m[35] += -1.5*kb/ODE_F.dt;
2171 }
2172 }
2173
2174 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2175 }
2176
2177
2178
2179
2180 void SMCZone::J3E_ddm_ombc_segment(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
2181 {
2182 Mat6 A;
2183 PetscInt index[6];
2184 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2185 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
2186 int z = zone_index;
2187 PetscScalar kb = mt->kb;
2188 PetscScalar e = mt->e;
2189 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2190 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2191 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2192 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2193 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2194 PetscScalar area = pcell->area;
2195 PetscScalar d_grad_T_dTi = 0;
2196 index[0] = zofs[z]+6*i+0;
2197 index[1] = zofs[z]+6*i+1;
2198 index[2] = zofs[z]+6*i+2;
2199 index[3] = zofs[z]+6*i+3;
2200 index[4] = zofs[z]+6*i+4;
2201 index[5] = zofs[z]+6*i+5;
2202
2203 for(int j=0;j<pcell->nb_num;j++)
2204 {
2205 int nb = pcell->nb_array[j];
2206 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2207 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
2208 //-------------------------------------
2209 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j]/area;
2210 //-------------------------------------
2211 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]/area; //dfun(3)/dT(r)
2212 if(j==0||j==pcell->nb_num-1)
2213 {
2214 PetscScalar h = pbc->heat_transfer;
2215 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
2216 d_grad_T_dTj += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
2217 }
2218 MatSetValue(*jac,zofs[z]+6*i+3,zofs[z]+6*nb+3,d_grad_T_dTj,INSERT_VALUES);
2219 }
2220
2221 Set_Mat6_I(A);
2222 A.m[21] = d_grad_T_dTi; //dfun(3)/dT(i)
2223
2224 A.m[25] = -Ti;
2225 A.m[27] = -ni;
2226
2227 A.m[32] = -Ti;
2228 A.m[33] = -pi;
2229 if(ODE_F.TimeDependent==true)
2230 {
2231 //second order
2232 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2233 {
2234 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2235 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2236 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2237 }
2238 else //first order
2239 {
2240 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2241 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2242 }
2243 }
2244 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2245
2246 int om_equ;
2247 int equ_num = 6;
2248 int size = pzone->davcell.size();
2249
2250 for(int j=0;j<electrode.size();j++)
2251 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
2252
2253 MatSetValue(*jac,zofs[zone_index]+6*i+0,zofs[zone_index]+equ_num*size+om_equ,-1,INSERT_VALUES);
2254 }
2255
2256
2257 void SMCZone::J3E_ddm_ombc_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc,
2258 ElZone *pz, int n)
2259 {
2260 Mat6 A;
2261 PetscInt index[6];
2262 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2263 OhmicBC *pbc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pcell->bc_index-1));
2264 int z = zone_index;
2265 PetscScalar kb = mt->kb;
2266 PetscScalar e = mt->e;
2267 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2268 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2269 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2270 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2271 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2272 PetscScalar area = pcell->area;
2273 PetscScalar d_grad_T_dTi = 0;
2274 index[0] = zofs[z]+6*i+0;
2275 index[1] = zofs[z]+6*i+1;
2276 index[2] = zofs[z]+6*i+2;
2277 index[3] = zofs[z]+6*i+3;
2278 index[4] = zofs[z]+6*i+4;
2279 index[5] = zofs[z]+6*i+5;
2280
2281 for(int j=0;j<pcell->nb_num;j++)
2282 {
2283 int nb = pcell->nb_array[j];
2284 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2285 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
2286 //-------------------------------------
2287 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j];
2288 //-------------------------------------
2289 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]; //dfun(3)/dT(r)
2290 MatSetValue(*jac,zofs[z]+6*i+3,zofs[z]+6*nb+3,d_grad_T_dTj,INSERT_VALUES);
2291 }
2292 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2293 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
2294 for(int j=0;j<ncell->nb_num;j++)
2295 {
2296 int nb = ncell->nb_array[j];
2297 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
2298 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
2299 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
2300 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
2301 MatSetValue(*jac,zofs[z]+6*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
2302 }
2303
2304 Set_Mat6_I(A);
2305 A.m[21] = d_grad_T_dTi; //dfun(3)/dT(i)
2306
2307 A.m[25] = -Ti;
2308 A.m[27] = -ni;
2309
2310 A.m[32] = -Ti;
2311 A.m[33] = -pi;
2312
2313 if(ODE_F.TimeDependent==true)
2314 {
2315 //second order
2316 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2317 {
2318 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2319 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2320 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2321 A.m[21] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
2322 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
2323 }
2324 else //first order
2325 {
2326 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2327 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2328 A.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
2329 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
2330 }
2331 }
2332
2333 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2334
2335 int om_equ;
2336 int equ_num = 6;
2337 int size = pzone->davcell.size();
2338
2339 for(int j=0;j<electrode.size();j++)
2340 if(electrode[j]==pcell->bc_index-1) { om_equ=j; break; }
2341
2342 MatSetValue(*jac,zofs[zone_index]+6*i+0,zofs[zone_index]+equ_num*size+om_equ,-1,INSERT_VALUES);
2343 }
2344
2345
2346 void SMCZone::J3E_ddm_stkbc_segment(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
2347 {
2348 Mat6 A;
2349 PetscInt index[6],col[6];
2350 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2351 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
2352 int z = zone_index;
2353 int equ_num = 6;
2354 int size = pzone->davcell.size();
2355
2356 PetscScalar kb = mt->kb;
2357 PetscScalar e = mt->e;
2358 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
2359 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
2360 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2361 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2362 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2363
2364 PetscScalar area = pcell->area;
2365 PetscScalar d_grad_T_dTi = 0;
2366 PetscScalar d_Fn_dni = 0;
2367 PetscScalar d_Fn_dTi = 0;
2368 PetscScalar d_Fp_dpi = 0;
2369 PetscScalar d_Fp_dTi = 0;
2370 //--------------------------------
2371 index[0] = zofs[z]+6*i+0;
2372 index[1] = zofs[z]+6*i+1;
2373 index[2] = zofs[z]+6*i+2;
2374 index[3] = zofs[z]+6*i+3;
2375 index[4] = zofs[z]+6*i+4;
2376 index[5] = zofs[z]+6*i+5;
2377 //--------------------------------
2378 for(int j=0;j<pcell->nb_num;j++)
2379 {
2380 int nb = pcell->nb_array[j];
2381 PetscScalar nj = x[zofs[z]+6*nb+1]; //electron density of nb node
2382 PetscScalar pj = x[zofs[z]+6*nb+2]; //hole density of nb node
2383 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2384 col[0] = zofs[z]+6*nb+0;
2385 col[1] = zofs[z]+6*nb+1;
2386 col[2] = zofs[z]+6*nb+2;
2387 col[3] = zofs[z]+6*nb+3;
2388 col[4] = zofs[z]+6*nb+4;
2389 col[5] = zofs[z]+6*nb+5;
2390 MatGetValues(*jtmp,6,index,6,col,A.m);
2391 A.m[0] = 0.0;
2392 A=A/area;
2393 if(j==0||j==pcell->nb_num-1)
2394 {
2395 d_Fn_dni += mt->band->pdSchottyJsn_pdn (ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2396 d_Fn_dTi += mt->band->pdSchottyJsn_pdTl(ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2397 d_Fp_dpi += mt->band->pdSchottyJsp_pdp (pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2398 d_Fp_dTi += mt->band->pdSchottyJsp_pdTl(pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2399 PetscScalar h = pbc->heat_transfer;
2400 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
2401 A.m[21] += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
2402 }
2403 A.m[24] = 0.0;
2404 A.m[25] = 0.0;
2405 A.m[26] = 0.0;
2406 A.m[27] = 0.0;
2407 A.m[28] = 0.0;
2408 A.m[29] = 0.0;
2409
2410 A.m[30] = 0.0;
2411 A.m[31] = 0.0;
2412 A.m[32] = 0.0;
2413 A.m[33] = 0.0;
2414 A.m[34] = 0.0;
2415 A.m[35] = 0.0;
2416 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2417 }
2418
2419 MatGetValues(*jtmp,6,index,6,index,A.m);
2420 A=A/area;
2421
2422 A.m[0] = 1.0;
2423
2424 A.m[7] += d_Fn_dni; //dfun(1)/dn(i)
2425 A.m[9] += d_Fn_dTi;
2426
2427 A.m[14]+= -d_Fp_dpi; //dfun(2)/dp(i)
2428 A.m[15]+= -d_Fp_dTi;
2429
2430 A.m[21]+= d_grad_T_dTi; //dfun(3)/dT(i)
2431
2432 A.m[24] = 0;
2433 A.m[25] = -Ti;
2434 A.m[26] = 0;
2435 A.m[27] = -ni;
2436 A.m[28] = 1.0;
2437 A.m[29] = 0;
2438
2439 A.m[30] = 0.0;
2440 A.m[31] = 0.0;
2441 A.m[32] = -Ti;
2442 A.m[33] = -pi;
2443 A.m[34] = 0.0;
2444 A.m[35] = 1.0;
2445 if(ODE_F.TimeDependent==true)
2446 {
2447 //second order
2448 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2449 {
2450 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2451 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2452 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2453 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2454 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2455 }
2456 else //first order
2457 {
2458 A.m[7] += -1/ODE_F.dt;
2459 A.m[14] += -1/ODE_F.dt;
2460 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2461 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2462 }
2463 }
2464 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2465
2466 int stk_equ;
2467 for(int j=0;j<electrode.size();j++)
2468 if(electrode[j]==pcell->bc_index-1)
2469 {
2470 stk_equ=j;
2471 break;
2472 }
2473 MatSetValue(*jac,zofs[zone_index]+6*i,zofs[zone_index]+equ_num*size+stk_equ,-1.0,INSERT_VALUES);
2474 }
2475
2476
2477 void SMCZone::J3E_ddm_stkbc_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc,
2478 ElZone *pz, int n)
2479 {
2480 Mat6 A;
2481 PetscInt index[6],col[6];
2482 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2483 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
2484 int z = zone_index;
2485 int equ_num = 6;
2486 int size = pzone->davcell.size();
2487 PetscScalar area = pcell->area;
2488 Vec6 scale;
2489 Set_Vec6(scale,area,area,area,1.0,area,area);
2490
2491 PetscScalar kb = mt->kb;
2492 PetscScalar e = mt->e;
2493 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
2494 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
2495 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2496 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2497 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2498
2499 PetscScalar d_grad_T_dTi = 0;
2500 PetscScalar d_Fn_dni = 0;
2501 PetscScalar d_Fn_dTi = 0;
2502 PetscScalar d_Fp_dpi = 0;
2503 PetscScalar d_Fp_dTi = 0;
2504 //--------------------------------
2505 index[0] = zofs[z]+6*i+0;
2506 index[1] = zofs[z]+6*i+1;
2507 index[2] = zofs[z]+6*i+2;
2508 index[3] = zofs[z]+6*i+3;
2509 index[4] = zofs[z]+6*i+4;
2510 index[5] = zofs[z]+6*i+5;
2511 //--------------------------------
2512 for(int j=0;j<pcell->nb_num;j++)
2513 {
2514 int nb = pcell->nb_array[j];
2515 PetscScalar nj = x[zofs[z]+6*nb+1]; //electron density of nb node
2516 PetscScalar pj = x[zofs[z]+6*nb+2]; //hole density of nb node
2517 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2518 col[0] = zofs[z]+6*nb+0;
2519 col[1] = zofs[z]+6*nb+1;
2520 col[2] = zofs[z]+6*nb+2;
2521 col[3] = zofs[z]+6*nb+3;
2522 col[4] = zofs[z]+6*nb+4;
2523 col[5] = zofs[z]+6*nb+5;
2524 MatGetValues(*jtmp,6,index,6,col,A.m);
2525 A.m[0] = 0.0;
2526 A=A/scale;
2527 if(j==0||j==pcell->nb_num-1)
2528 {
2529 d_Fn_dni += mt->band->pdSchottyJsn_pdn (ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2530 d_Fn_dTi += mt->band->pdSchottyJsn_pdTl(ni,Ti,VB-deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2531 d_Fp_dpi += mt->band->pdSchottyJsp_pdp (pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2532 d_Fp_dTi += mt->band->pdSchottyJsp_pdTl(pi,Ti,VB+deltaVB)*0.5*pcell->ilen[j]/pcell->area;
2533 }
2534 A.m[24] = 0.0;
2535 A.m[25] = 0.0;
2536 A.m[26] = 0.0;
2537 A.m[27] = 0.0;
2538 A.m[28] = 0.0;
2539 A.m[29] = 0.0;
2540
2541 A.m[30] = 0.0;
2542 A.m[31] = 0.0;
2543 A.m[32] = 0.0;
2544 A.m[33] = 0.0;
2545 A.m[34] = 0.0;
2546 A.m[35] = 0.0;
2547 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2548 }
2549
2550 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2551 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
2552 for(int j=0;j<ncell->nb_num;j++)
2553 {
2554 int nb = ncell->nb_array[j];
2555 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
2556 PetscScalar kapa = pz->mt->thermal->HeatConduction(fs[i].T);//(0.5*(Ti+Tj_n));
2557 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
2558 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
2559 MatSetValue(*jac,zofs[z]+6*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
2560 }
2561
2562 MatGetValues(*jtmp,6,index,6,index,A.m);
2563 A=A/scale;
2564
2565 A.m[0] = 1.0;
2566
2567 A.m[7] += d_Fn_dni; //dfun(1)/dn(i)
2568 A.m[9] += d_Fn_dTi;
2569
2570 A.m[14]+= -d_Fp_dpi; //dfun(2)/dp(i)
2571 A.m[15]+= -d_Fp_dTi;
2572
2573 A.m[21]+= d_grad_T_dTi; //dfun(3)/dT(i)
2574
2575 A.m[24] = 0;
2576 A.m[25] = -Ti;
2577 A.m[26] = 0;
2578 A.m[27] = -ni;
2579 A.m[28] = 1.0;
2580 A.m[29] = 0;
2581
2582 A.m[30] = 0.0;
2583 A.m[31] = 0.0;
2584 A.m[32] = -Ti;
2585 A.m[33] = -pi;
2586 A.m[34] = 0.0;
2587 A.m[35] = 1.0;
2588
2589 if(ODE_F.TimeDependent==true)
2590 {
2591 //second order
2592 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2593 {
2594 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2595 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2596 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2597 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2598 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2599 A.m[21] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
2600 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
2601 }
2602 else //first order
2603 {
2604 A.m[7] += -1/ODE_F.dt;
2605 A.m[14] += -1/ODE_F.dt;
2606 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2607 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2608 A.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
2609 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
2610 }
2611 }
2612
2613 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2614
2615 int stk_equ;
2616 for(int j=0;j<electrode.size();j++)
2617 if(electrode[j]==pcell->bc_index-1)
2618 {
2619 stk_equ=j;
2620 break;
2621 }
2622 MatSetValue(*jac,zofs[zone_index]+6*i,zofs[zone_index]+equ_num*size+stk_equ,-1.0,INSERT_VALUES);
2623 }
2624
2625
2626 void SMCZone::J3E_ddm_insulator_gate(int i,PetscScalar *x, Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
2627 {
2628 Mat6 A;
2629 PetscInt index[6],col[6];
2630 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2631 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
2632 int z = zone_index;
2633 int equ_num = 6;
2634 int size = pzone->davcell.size();
2635 PetscScalar kb = mt->kb;
2636 PetscScalar e = mt->e;
2637 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2638 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2639 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2640 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2641 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2642 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2643
2644 PetscScalar area = pcell->area;
2645 PetscScalar d_grad_P_dVi = 0;
2646 PetscScalar d_grad_P_dVapp = 0;
2647 PetscScalar d_grad_T_dTi = 0;
2648 int ins_equ;
2649 for(int j=0;j<electrode.size();j++)
2650 if(electrode[j]==pcell->bc_index-1)
2651 {ins_equ=j;break;}
2652 //--------------------------------
2653 index[0] = zofs[z]+6*i+0;
2654 index[1] = zofs[z]+6*i+1;
2655 index[2] = zofs[z]+6*i+2;
2656 index[3] = zofs[z]+6*i+3;
2657 index[4] = zofs[z]+6*i+4;
2658 index[5] = zofs[z]+6*i+5;
2659 //--------------------------------
2660 for(int j=0;j<pcell->nb_num;j++)
2661 {
2662 int nb = pcell->nb_array[j];
2663 PetscScalar Vj = x[zofs[z]+6*nb+0]; //potential of nb node
2664 PetscScalar Tj = x[zofs[z]+6*nb+3]; //lattice temperature of nb node
2665 col[0] = zofs[z]+6*nb+0;
2666 col[1] = zofs[z]+6*nb+1;
2667 col[2] = zofs[z]+6*nb+2;
2668 col[3] = zofs[z]+6*nb+3;
2669 col[4] = zofs[z]+6*nb+4;
2670 col[5] = zofs[z]+6*nb+5;
2671
2672 MatGetValues(*jtmp,6,index,6,col,A.m);
2673 A=A/area;
2674 if(j==0||j==pcell->nb_num-1)
2675 {
2676 PetscScalar Thick = pbc->Thick;
2677 PetscScalar eps_ox = mt->eps0*pbc->eps;
2678 PetscScalar eps = aux[i].eps;
2679 PetscScalar s=eps_ox/eps/Thick;
2680 d_grad_P_dVi += -0.5*pcell->ilen[j]*0.25*s*3/area;
2681 d_grad_P_dVapp += 0.5*pcell->ilen[j]*s/area;
2682 A.m[0]+= -0.5*pcell->ilen[j]*0.25*s/area;
2683 PetscScalar h = pbc->heat_transfer;
2684 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/area;
2685 A.m[21] += -0.5*pcell->ilen[j]*0.25*h/area;
2686 }
2687 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2688
2689 }
2690 //fun(0) is the poisson's equation
2691 //fun(1) is the continuous equation of electron
2692 //fun(2) is the continuous equation of hole
2693
2694 MatGetValues(*jtmp,6,index,6,index,A.m);
2695 A=A/area;
2696 A.m[0] += d_grad_P_dVi;
2697 A.m[1] += -e/aux[i].eps; //df(phi)/dn(i)
2698 A.m[2] += e/aux[i].eps; //df(phi)/dp(i)
2699
2700 A.m[21]+= d_grad_T_dTi;
2701
2702 if(ODE_F.TimeDependent==true)
2703 {
2704 //second order
2705 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2706 {
2707 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2708 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2709 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2710 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2711 A.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2712 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2713 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2714 }
2715 else //first order
2716 {
2717 A.m[7] += -1/ODE_F.dt;
2718 A.m[14] += -1/ODE_F.dt;
2719 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2720 A.m[21] += -aux[i].density*HeatCapacity/ODE_F.dt;
2721 A.m[28] += -1.5*kb/ODE_F.dt;
2722 A.m[35] += -1.5*kb/ODE_F.dt;
2723 }
2724 }
2725 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2726 MatSetValue(*jac,zofs[z]+6*i,zofs[z]+equ_num*size+ins_equ,d_grad_P_dVapp,INSERT_VALUES);
2727
2728 }
2729
2730
2731 void SMCZone::J3E_ddm_interface(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
2732 ISZone *pz, int n)
2733 {
2734 Mat6 A;
2735 PetscInt index[6],col[6];
2736 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2737 InsulatorInterfaceBC *pbc = dynamic_cast<InsulatorInterfaceBC * >(bc.Get_pointer(pcell->bc_index-1));
2738 int z = zone_index;
2739 PetscScalar kb = mt->kb;
2740 PetscScalar e = mt->e;
2741 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2742 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2743 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2744 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2745 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2746 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2747
2748 PetscScalar area = pcell->area;
2749 PetscScalar d_grad_P_dVi = 0;
2750 PetscScalar d_grad_T_dTi = 0;
2751 Vec6 scale;
2752 Set_Vec6(scale,1.0,area,area,1.0,area,area);
2753 //--------------------------------
2754 index[0] = zofs[z]+6*i+0;
2755 index[1] = zofs[z]+6*i+1;
2756 index[2] = zofs[z]+6*i+2;
2757 index[3] = zofs[z]+6*i+3;
2758 index[4] = zofs[z]+6*i+4;
2759 index[5] = zofs[z]+6*i+5;
2760 //--------------------------------
2761 for(int j=0;j<pcell->nb_num;j++)
2762 {
2763 int nb = pcell->nb_array[j];
2764
2765 PetscScalar Vj = x[zofs[z]+6*nb+0]; //potential of nb node
2766 PetscScalar Tj = x[zofs[z]+6*nb+3];
2767 PetscScalar kapa = mt->thermal->HeatConduction(0.5*(Ti+Tj));
2768 d_grad_P_dVi += -aux[i].eps*pcell->elen[j]/pcell->ilen[j];
2769 col[0] = zofs[z]+6*nb+0;
2770 col[1] = zofs[z]+6*nb+1;
2771 col[2] = zofs[z]+6*nb+2;
2772 col[3] = zofs[z]+6*nb+3;
2773 col[4] = zofs[z]+6*nb+4;
2774 col[5] = zofs[z]+6*nb+5;
2775 MatGetValues(*jtmp,6,index,6,col,A.m);
2776 A=A/scale;
2777 A.m[0] = aux[i].eps*pcell->elen[j]/pcell->ilen[j]; //dfun(0)/dP(r)
2778 A.m[1] = 0; //dfun(0)/dn(r)
2779 A.m[2] = 0; //dfun(0)/dp(r)
2780 A.m[3] = 0; //dfun(0)/dT(r)
2781 A.m[4] = 0;
2782 A.m[5] = 0;
2783 MatSetValues(*jac,6,index,6,col,A.m,INSERT_VALUES);
2784 }
2785
2786 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2787 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
2788 for(int j=0;j<ncell->nb_num;j++)
2789 {
2790 int nb = ncell->nb_array[j];
2791 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
2792 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
2793 d_grad_P_dVi += -pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
2794 PetscScalar d_grad_P_dVj = pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
2795 MatSetValue(*jac,zofs[z]+6*i+0,zofs[pz->pzone->zone_index]+2*nb+0,d_grad_P_dVj,INSERT_VALUES);
2796 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
2797 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
2798 MatSetValue(*jac,zofs[z]+6*i+3,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
2799 }
2800
2801 //fun(0) is the poisson's equation
2802 //fun(1) is the continuous equation of electron
2803 //fun(2) is the continuous equation of hole
2804
2805 MatGetValues(*jtmp,6,index,6,index,A.m);
2806 A=A/scale;
2807 A.m[0] = d_grad_P_dVi; //dfun(0)/dP(i)
2808 A.m[1] = -e*pcell->area; //dfun(0)/dn(i)
2809 A.m[2] = e*pcell->area; //dfun(0)/dp(i)
2810 A.m[3] = 0.0; //dfun(0)/dT(i)
2811 A.m[4] = 0.0;
2812 A.m[5] = 0.0;
2813
2814 A.m[21]+= d_grad_T_dTi;
2815
2816
2817 if(ODE_F.TimeDependent==true)
2818 {
2819 //second order
2820 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2821 {
2822 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2823 A.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2824 A.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2825 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2826 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2827 A.m[21] += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
2828 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
2829 A.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2830 A.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2831 }
2832 else //first order
2833 {
2834 A.m[7] += -1/ODE_F.dt;
2835 A.m[14] += -1/ODE_F.dt;
2836 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2837 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
2838 A.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
2839 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
2840 A.m[28] += -1.5*kb/ODE_F.dt;
2841 A.m[35] += -1.5*kb/ODE_F.dt;
2842 }
2843 }
2844 MatSetValues(*jac,6,index,6,index,A.m,INSERT_VALUES);
2845 }
2846
2847
2848
2849 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)
2850 {
2851 Mat6 A1,A2,AJ;
2852 PetscInt index[6],indexn[6],col[6];
2853 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
2854 int z = zone_index;
2855 PetscScalar kb = mt->kb;
2856 PetscScalar e = mt->e;
2857 PetscScalar Vi = x[zofs[z]+6*i+0]; //potential of node i
2858 PetscScalar ni = x[zofs[z]+6*i+1]; //electron density of node i
2859 PetscScalar pi = x[zofs[z]+6*i+2]; //hole density of node i
2860 PetscScalar Ti = x[zofs[z]+6*i+3]; //lattice temperature of node i
2861 PetscScalar Tn = x[zofs[z]+6*i+4]/ni;
2862 PetscScalar Tp = x[zofs[z]+6*i+5]/pi;
2863 //--------------------------------
2864 index[0] = zofs[z]+6*i+0;
2865 index[1] = zofs[z]+6*i+1;
2866 index[2] = zofs[z]+6*i+2;
2867 index[3] = zofs[z]+6*i+3;
2868 index[4] = zofs[z]+6*i+4;
2869 index[5] = zofs[z]+6*i+5;
2870 //--------------------------------
2871 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
2872 indexn[0] = zofs[pz->zone_index]+6*n+0;
2873 indexn[1] = zofs[pz->zone_index]+6*n+1;
2874 indexn[2] = zofs[pz->zone_index]+6*n+2;
2875 indexn[3] = zofs[pz->zone_index]+6*n+3;
2876 indexn[4] = zofs[pz->zone_index]+6*n+4;
2877 indexn[5] = zofs[pz->zone_index]+6*n+5;
2878 //--------------------------------
2879 PetscScalar area = pcell->area + ncell->area;
2880 //--------------------------------
2881 if(zone_index > pz->pzone->zone_index)
2882 {
2883 Set_Mat6_I(AJ);
2884 MatSetValues(*jac,6,index,6,index,AJ.m,INSERT_VALUES);
2885 AJ = -1.0*AJ;
2886 MatSetValues(*jac,6,index,6,indexn,AJ.m,INSERT_VALUES);
2887 return;
2888 }
2889 //--------------------------------
2890 for(int j=0;j<pcell->nb_num;j++)
2891 {
2892 int nb = pcell->nb_array[j];
2893 col[0] = zofs[z]+6*nb+0;
2894 col[1] = zofs[z]+6*nb+1;
2895 col[2] = zofs[z]+6*nb+2;
2896 col[3] = zofs[z]+6*nb+3;
2897 col[4] = zofs[z]+6*nb+4;
2898 col[5] = zofs[z]+6*nb+5;
2899 MatGetValues(*jtmp,6,index,6,col,A1.m);
2900 //-------------------------------------
2901 A1=A1/area;
2902 MatSetValues(*jac,6,index,6,col,A1.m,INSERT_VALUES);
2903 }
2904
2905 //process another half cell in dornor zone
2906 for(int j=0;j<ncell->nb_num;j++)
2907 {
2908 int nb = ncell->nb_array[j];
2909 col[0] = zofs[pz->zone_index]+6*nb+0;
2910 col[1] = zofs[pz->zone_index]+6*nb+1;
2911 col[2] = zofs[pz->zone_index]+6*nb+2;
2912 col[3] = zofs[pz->zone_index]+6*nb+3;
2913 col[4] = zofs[pz->zone_index]+6*nb+4;
2914 col[5] = zofs[pz->zone_index]+6*nb+5;
2915 MatGetValues(*jtmp,6,indexn,6,col,A2.m);
2916 A2=A2/area;
2917 MatSetValues(*jac,6,index,6,col,A2.m,INSERT_VALUES);
2918 }
2919
2920 //fun(0) is the poisson's equation
2921 //fun(1) is the continuous equation of electron
2922 //fun(2) is the continuous equation of hole
2923 MatGetValues(*jtmp,6,index,6,index,A1.m);
2924 MatGetValues(*jtmp,6,indexn,6,indexn,A2.m);
2925 AJ=(A1+A2)/area;
2926
2927 AJ.m[1] += -e/aux[i].eps; //dfun(0)/dn(i)
2928 AJ.m[2] += e/aux[i].eps; //dfun(0)/dp(i)
2929
2930 if(ODE_F.TimeDependent==true)
2931 {
2932 //second order
2933 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
2934 {
2935 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
2936 AJ.m[7] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2937 AJ.m[14] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2938 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
2939 AJ.m[21] += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2940 AJ.m[28] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2941 AJ.m[35] += -1.5*kb*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
2942 }
2943 else //first order
2944 {
2945 AJ.m[7] += -1/ODE_F.dt;
2946 AJ.m[14] += -1/ODE_F.dt;
2947 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
2948 AJ.m[21] += -aux[i].density*HeatCapacity1/ODE_F.dt;
2949 AJ.m[28] += -1.5*kb/ODE_F.dt;
2950 AJ.m[35] += -1.5*kb/ODE_F.dt;
2951 }
2952 }
2953 MatSetValues(*jac,6,index,6,index,AJ.m,INSERT_VALUES);
2954 }
2955
2956
2957 void SMCZone::J3E_om_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
2958 {
2959 int equ_num = 6;
2960 int size = pzone->davcell.size();
2961 int bc_index = electrode[i];
2962 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
2963 PetscScalar A1[6],A2[6],J[6];
2964 PetscInt index[6],col[6];
2965 int matrix_row = zofs[zone_index]+equ_num*size+i;
2966 int connect_index;
2967 if(pbc->inner_connect!=-1)
2968 {
2969 int connect_to = pbc->inner_connect;
2970 int connect_zone = pbc->connect_zone;
2971 int connect_elec=-1;
2972 SMCZone * pz = dynamic_cast <SMCZone * > (pbc->pzonedata);
2973 for(int j=0;j<pz->electrode.size();j++)
2974 {
2975 if(bc[pz->electrode[j]].BCType==OhmicContact)
2976 {
2977 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pz->electrode[j]));
2978 if(om_bc->inner_connect==bc_index) connect_elec=j;
2979 }
2980 }
2981 connect_index=zofs[connect_zone]+equ_num*pz->pzone->davcell.size()+connect_elec;
2982 }
2983 PetscScalar R = pbc->R;
2984 PetscScalar C = pbc->C;
2985 PetscScalar L = pbc->L;
2986 PetscScalar kb = mt->kb;
2987 PetscScalar e = mt->e;
2988
2989 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
2990 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
2991
2992 for(int n=0;n<bc[bc_index].psegment->node_array.size();n++)
2993 {
2994 int node=bc[bc_index].psegment->node_array[n];
2995 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
2996 PetscScalar dJdisp_dVi=0,dJdisp_dVr=0;
2997 index[0] = zofs[zone_index]+6*node+0;
2998 index[1] = zofs[zone_index]+6*node+1;
2999 index[2] = zofs[zone_index]+6*node+2;
3000 index[3] = zofs[zone_index]+6*node+3;
3001 index[4] = zofs[zone_index]+6*node+4;
3002 index[5] = zofs[zone_index]+6*node+5;
3003 for(int j=0;j<pcell->nb_num;j++)
3004 {
3005 int nb = pcell->nb_array[j];
3006 col[0] = zofs[zone_index]+6*nb+0;
3007 col[1] = zofs[zone_index]+6*nb+1;
3008 col[2] = zofs[zone_index]+6*nb+2;
3009 col[3] = zofs[zone_index]+6*nb+3;
3010 col[4] = zofs[zone_index]+6*nb+4;
3011 col[5] = zofs[zone_index]+6*nb+5;
3012
3013 MatGetValues(*jtmp,1,&index[1],6,col,A1);
3014 MatGetValues(*jtmp,1,&index[2],6,col,A2);
3015
3016 //for displacement current
3017 dJdisp_dVi += aux[node].eps/pcell->ilen[j]/ODE_F.dt*pcell->elen[j];
3018 dJdisp_dVr = -aux[node].eps/pcell->ilen[j]/ODE_F.dt*pcell->elen[j];
3019 if(pbc->inner_connect!=-1)
3020 {
3021 int connect_to = pbc->inner_connect;
3022 if(bc_index < connect_to)
3023 {
3024 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3025 J[1]=(A1[1]-A2[1])*DeviceDepth;
3026 J[2]=(A1[2]-A2[2])*DeviceDepth;
3027 J[3]=(A1[3]-A2[3])*DeviceDepth;
3028 J[4]=(A1[4]-A2[4])*DeviceDepth;
3029 J[5]=(A1[5]-A2[5])*DeviceDepth;
3030 MatSetValues(*jac,1,&connect_index,6,col,J,ADD_VALUES);
3031 }
3032 else
3033 {
3034 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3035 J[1]=(A1[1]-A2[1])*DeviceDepth;
3036 J[2]=(A1[2]-A2[2])*DeviceDepth;
3037 J[3]=(A1[3]-A2[3])*DeviceDepth;
3038 J[4]=(A1[4]-A2[4])*DeviceDepth;
3039 J[5]=(A1[5]-A2[5])*DeviceDepth;
3040 MatSetValues(*jac,1,&matrix_row,6,col,J,ADD_VALUES);
3041 }
3042 }
3043 else if(pbc->electrode_type==VoltageBC)
3044 {
3045 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth*(L/ODE_F.dt+R);
3046 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
3047 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
3048 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
3049 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
3050 J[5]=(A1[5]-A2[5])*DeviceDepth*(L/ODE_F.dt+R);
3051 MatSetValues(*jac,1,&matrix_row,6,col,J,ADD_VALUES);
3052 }
3053 else if(pbc->electrode_type==CurrentBC)
3054 {
3055 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3056 J[1]=(A1[1]-A2[1])*DeviceDepth;
3057 J[2]=(A1[2]-A2[2])*DeviceDepth;
3058 J[3]=(A1[3]-A2[3])*DeviceDepth;
3059 J[4]=(A1[4]-A2[4])*DeviceDepth;
3060 J[5]=(A1[5]-A2[5])*DeviceDepth;
3061 MatSetValues(*jac,1,&matrix_row,6,col,J,ADD_VALUES);
3062 }
3063 }
3064 MatGetValues(*jtmp,1,&index[1],6,index,A1);
3065 MatGetValues(*jtmp,1,&index[2],6,index,A2);
3066 if(pbc->inner_connect!=-1)
3067 {
3068 int connect_to = pbc->inner_connect;
3069 if(bc_index < connect_to)
3070 {
3071 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3072 J[1]=(A1[1]-A2[1])*DeviceDepth;
3073 J[2]=(A1[2]-A2[2])*DeviceDepth;
3074 J[3]=(A1[3]-A2[3])*DeviceDepth;
3075 J[4]=(A1[4]-A2[4])*DeviceDepth;
3076 J[5]=(A1[5]-A2[5])*DeviceDepth;
3077 MatSetValues(*jac,1,&connect_index,6,index,J,ADD_VALUES);
3078 }
3079 else
3080 {
3081 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
3082 J[1]=(A1[1]-A2[1])*DeviceDepth;
3083 J[2]=(A1[2]-A2[2])*DeviceDepth;
3084 J[3]=(A1[3]-A2[3])*DeviceDepth;
3085 J[4]=(A1[4]-A2[4])*DeviceDepth;
3086 J[5]=(A1[5]-A2[5])*DeviceDepth;
3087 MatSetValues(*jac,1,&matrix_row,6,index,J,ADD_VALUES);
3088 }
3089 }
3090 else if(pbc->electrode_type==VoltageBC)
3091 {
3092 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth*(L/ODE_F.dt+R);
3093 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
3094 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
3095 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
3096 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
3097 J[5]=(A1[5]-A2[5])*DeviceDepth*(L/ODE_F.dt+R);
3098 MatSetValues(*jac,1,&matrix_row,6,index,J,ADD_VALUES);
3099 }
3100 else if(pbc->electrode_type==CurrentBC)
3101 {
3102 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth;
3103 J[1]=(A1[1]-A2[1])*DeviceDepth;
3104 J[2]=(A1[2]-A2[2])*DeviceDepth;
3105 J[3]=(A1[3]-A2[3])*DeviceDepth;
3106 J[4]=(A1[4]-A2[4])*DeviceDepth;
3107 J[5]=(A1[5]-A2[5])*DeviceDepth;
3108 MatSetValues(*jac,1,&matrix_row,6,index,J,ADD_VALUES);
3109 }
3110 }
3111
3112 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3113 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3114 if(pbc->inner_connect!=-1)
3115 {
3116 int connect_to = pbc->inner_connect;
3117 if(bc_index < connect_to)
3118 {
3119 MatSetValue(*jac,matrix_row,matrix_row,1.0,ADD_VALUES);
3120 }
3121 else
3122 {
3123 MatSetValue(*jac,connect_index,matrix_row,-1.0,ADD_VALUES);
3124 }
3125 }
3126 else if(pbc->electrode_type==VoltageBC)
3127 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
3128
3129 }
3130
3131
3132 void SMCZone::J3E_stk_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
3133 {
3134 int equ_num = 6;
3135 int size = pzone->davcell.size();
3136 int bc_index = electrode[i];
3137 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
3138 int matrix_row = zofs[zone_index]+equ_num*size+i;
3139 PetscScalar R = pbc->R;
3140 PetscScalar C = pbc->C;
3141 PetscScalar L = pbc->L;
3142 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
3143 {
3144 int node = bc[bc_index].psegment->node_array[j];
3145 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
3146 PetscScalar VB = pbc->WorkFunction-aux[node].affinity;
3147 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
3148 PetscScalar ni = x[zofs[zone_index]+6*node+1]; //electron density of node i
3149 PetscScalar pi = x[zofs[zone_index]+6*node+2]; //hole density of node i
3150 PetscScalar dJ_dni = -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[0]
3151 -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
3152 PetscScalar dJ_dpi = -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[0]
3153 -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
3154
3155 //for displacement current
3156 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3157 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3158 for(int k=0;k<pcell->nb_num;k++)
3159 {
3160 int nb = pcell->nb_array[k];
3161 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3162 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3163 if(pbc->electrode_type==VoltageBC)
3164 {
3165 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+0,dI_dVi*(L/ODE_F.dt+R),ADD_VALUES);
3166 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*nb+0,dI_dVr*(L/ODE_F.dt+R),ADD_VALUES);
3167 }
3168 else if(pbc->electrode_type==CurrentBC)
3169 {
3170 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+0,dI_dVi,ADD_VALUES);
3171 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*nb+0,dI_dVr,ADD_VALUES);
3172 }
3173 }
3174 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3175 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3176
3177 if(pbc->electrode_type==VoltageBC)
3178 {
3179 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+1,DeviceDepth*dJ_dni*(L/ODE_F.dt+R),INSERT_VALUES);
3180 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+2,DeviceDepth*dJ_dpi*(L/ODE_F.dt+R),INSERT_VALUES);
3181 }
3182 else if(pbc->electrode_type==CurrentBC)
3183 {
3184 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+1,DeviceDepth*dJ_dni,INSERT_VALUES);
3185 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+2,DeviceDepth*dJ_dpi,INSERT_VALUES);
3186 }
3187 }
3188 if(pbc->electrode_type==VoltageBC)
3189 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
3190 }
3191
3192
3193 void SMCZone::J3E_ins_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
3194 {
3195 int equ_num = 6;
3196 int size = pzone->davcell.size();
3197 int bc_index = electrode[i];
3198 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
3199 int matrix_row = zofs[zone_index]+equ_num*size+i;
3200 PetscScalar R = pbc->R;
3201 PetscScalar C = pbc->C;
3202 PetscScalar L = pbc->L;
3203
3204 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3205 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3206 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
3207 {
3208 int node = bc[bc_index].psegment->node_array[j];
3209 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
3210
3211 //for displacement current
3212
3213 for(int k=0;k<pcell->nb_num;k++)
3214 {
3215 int nb = pcell->nb_array[k];
3216 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3217 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
3218 if(pbc->electrode_type==VoltageBC)
3219 {
3220 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*node+0,dI_dVi*(L/ODE_F.dt+R),ADD_VALUES);
3221 MatSetValue(*jac,matrix_row,zofs[zone_index]+6*nb+0,dI_dVr*(L/ODE_F.dt+R),ADD_VALUES);
3222 }
3223 }
3224
3225 }
3226 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
3227 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
3228
3229 if(pbc->electrode_type==VoltageBC)
3230 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
3231 }
3232
3233
3234 void SMCZone::F3E_efield_update(PetscScalar *x,vector<int> & zofs, DABC &bc, vector<BZoneData *>zonedata)
3235 {
3236 PetscScalar P_x=0,P_y=0,w=0;
3237 //calculate Ex Ey with least-squares gradient construction
3238 for(int i=0; i<pzone->davcell.size();i++)
3239 {
3240 P_x=0,P_y=0,w=0;
3241 VoronoiCell *pcell = pzone->davcell.GetPointer(i);
3242 for(int k=0;k<pcell->nb_num;k++)
3243 {
3244 const VoronoiCell *ncell = pzone->davcell.GetPointer(pcell->nb_array[k]);
3245 PetscScalar dx = ncell->x - pcell->x;
3246 PetscScalar dy = ncell->y - pcell->y;
3247 PetscScalar dP = x[zofs[zone_index]+6*pcell->nb_array[k]+0] - x[zofs[zone_index]+6*i+0];
3248 w=1.0/sqrt(dx*dx+dy*dy);
3249 P_x+=w*w*dP*dx;
3250 P_y+=w*w*dP*dy;
3251 }
3252 if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
3253 {
3254 InsulatorInterfaceBC *pbc;
3255 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
3256 int n_zone = pbc->pinterface->Find_neighbor_zone_index(zone_index);
3257 int n_node = pbc->pinterface->Find_neighbor_node_index(zone_index,i);
3258 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
3259 const VoronoiCell* dcell = pz->pzone->davcell.GetPointer(n_node);
3260 for(int k=1;k<dcell->nb_num-1;k++)
3261 {
3262 const VoronoiCell *ncell = pz->pzone->davcell.GetPointer(dcell->nb_array[k]);
3263 PetscScalar dx = ncell->x - dcell->x;
3264 PetscScalar dy = ncell->y - dcell->y;
3265 PetscScalar dP = x[zofs[n_zone]+2*dcell->nb_array[k]]- x[zofs[n_zone]+2*n_node];
3266 w=1.0/sqrt(dx*dx+dy*dy);
3267 P_x+=w*w*dP*dx;
3268 P_y+=w*w*dP*dy;
3269 }
3270 }
3271 aux[i].Ex = -(pcell->sc*P_x-pcell->sb*P_y);
3272 aux[i].Ey = -(pcell->sa*P_y-pcell->sb*P_x);
3273 }
3274 }
General-purpose Semiconductor Simulator