?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
more fix on Ec/Ev.
[gss-tcad.git] / src / solver / qddm1e / semiequ1q.cc
blob6f24182e8d8e56416265058bb3a471ab15901531
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: July 20, 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 "zonedata.h"
22 #include "jflux1q.h"
23 #include "vec5.h"
24 #define _FLUX1_
25 #define _HMOB_
26 #define _EXP_R_
27
28
29 void SMCZone::F1Q_Tri_qddm(Tri *ptri,PetscScalar *x,PetscScalar *f, vector<int> &zofs)
30 {
31 int A = ptri->node[0];
32 int B = ptri->node[1];
33 int C = ptri->node[2];
34 PetscScalar xa = pzone->danode[ptri->node[0]].x;
35 PetscScalar xb = pzone->danode[ptri->node[1]].x;
36 PetscScalar xc = pzone->danode[ptri->node[2]].x;
37 PetscScalar ya = pzone->danode[ptri->node[0]].y;
38 PetscScalar yb = pzone->danode[ptri->node[1]].y;
39 PetscScalar yc = pzone->danode[ptri->node[2]].y;
40 PetscScalar tri_area = ptri->area;
41 PetscScalar Sax = (xc-xb)/ptri->edge_len[0];
42 PetscScalar Say = (yc-yb)/ptri->edge_len[0];
43 PetscScalar Sbx = (xa-xc)/ptri->edge_len[1];
44 PetscScalar Sby = (ya-yc)/ptri->edge_len[1];
45 PetscScalar Scx = (xb-xa)/ptri->edge_len[2];
46 PetscScalar Scy = (yb-ya)/ptri->edge_len[2];
47 PetscScalar Wab = ptri->d[1]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
48 PetscScalar Wac = ptri->d[2]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[1])*sin(ptri->angle[1]));
49 PetscScalar Wbc = ptri->d[2]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[0])*sin(ptri->angle[0]));
50 PetscScalar Wba = ptri->d[0]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
51 PetscScalar Wca = ptri->d[0]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[1])*sin(ptri->angle[1]));
52 PetscScalar Wcb = ptri->d[1]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[0])*sin(ptri->angle[0]));
53 PetscScalar Ma = -cos(ptri->angle[0]);
54 PetscScalar Mb = -cos(ptri->angle[1]);
55 PetscScalar Mc = -cos(ptri->angle[2]);
56
57 PetscScalar e = mt->e;
58 PetscScalar hbar = mt->hbar;
59 PetscScalar T = (fs[A].T+fs[B].T+fs[C].T)/3.0;
60 PetscScalar kb = mt->kb;
61 PetscScalar Vt = mt->kb*T/e;
62 PetscScalar Eg = mt->band->Eg(T);
63 PetscScalar me = mt->band->EffecElecMass(T);
64 PetscScalar mh = mt->band->EffecHoleMass(T);
65 PetscScalar gamman = mt->band->Gamman();
66 PetscScalar gammap = mt->band->Gammap();
67 PetscScalar bn = gamman*hbar*hbar/(6*e*me);
68 PetscScalar bp = gammap*hbar*hbar/(6*e*mh);
69 PetscScalar S;
70
71 PetscScalar Va = x[zofs[zone_index]+5*A+0]; //potential of node A
72 PetscScalar na = x[zofs[zone_index]+5*A+1]; //electron density of node A
73 PetscScalar pa = x[zofs[zone_index]+5*A+2]; //hole density of node A
74 PetscScalar Eqca = x[zofs[zone_index]+5*A+3]; //quantum conduction band of node A
75 PetscScalar Eqva = x[zofs[zone_index]+5*A+4]; //quantum valence band of node A
76 mt->mapping(&pzone->danode[A],&aux[A],0);
77 PetscScalar Eca = e*(Eqca/e + kb*fs[A].T/e*(log(aux[A].Nc)-1.5*log(fs[A].T)));
78 PetscScalar Eva = e*(Eqva/e - kb*fs[A].T/e*(log(aux[A].Nv)-1.5*log(fs[A].T)));
79 PetscScalar Ra = mt->band->Recomb(pa,na,fs[A].T);
80 PetscScalar nia = mt->band->nie(T);
81 PetscScalar phina = Va - log(fabs(na)/nia)*Vt;
82 PetscScalar phipa = Va + log(fabs(pa)/nia)*Vt;
83 PetscScalar qca = -(-e*phina-Eqca)/(kb*T);
84 PetscScalar qva = (-e*phipa-Eqva)/(kb*T);
85
86 PetscScalar Vb = x[zofs[zone_index]+5*B+0]; //potential of node B
87 PetscScalar nb = x[zofs[zone_index]+5*B+1]; //electron density of node B
88 PetscScalar pb = x[zofs[zone_index]+5*B+2]; //hole density of node B
89 PetscScalar Eqcb = x[zofs[zone_index]+5*B+3]; //quantum conduction band of node B
90 PetscScalar Eqvb = x[zofs[zone_index]+5*B+4]; //quantum valence band of node B
91 mt->mapping(&pzone->danode[B],&aux[B],0);
92 PetscScalar Ecb = e*(Eqcb/e + kb*fs[B].T/e*(log(aux[B].Nc)-1.5*log(fs[B].T)));
93 PetscScalar Evb = e*(Eqvb/e - kb*fs[B].T/e*(log(aux[B].Nv)-1.5*log(fs[B].T)));
94 PetscScalar Rb = mt->band->Recomb(pb,nb,fs[B].T);
95 PetscScalar nib = mt->band->nie(T);
96 PetscScalar phinb = Vb - log(fabs(nb)/nib)*Vt;
97 PetscScalar phipb = Vb + log(fabs(pb)/nib)*Vt;
98 PetscScalar qcb = -(-e*phinb-Eqcb)/(kb*T);
99 PetscScalar qvb = (-e*phipb-Eqvb)/(kb*T);
100
101 PetscScalar Vc = x[zofs[zone_index]+5*C+0]; //potential of node C
102 PetscScalar nc = x[zofs[zone_index]+5*C+1]; //electron density of node C
103 PetscScalar pc = x[zofs[zone_index]+5*C+2]; //hole density of node C
104 PetscScalar Eqcc = x[zofs[zone_index]+5*C+3]; //quantum conduction band of node C
105 PetscScalar Eqvc = x[zofs[zone_index]+5*C+4]; //quantum valence band of node C
106 mt->mapping(&pzone->danode[C],&aux[C],0);
107 PetscScalar Ecc = e*(Eqcc/e + kb*fs[C].T/e*(log(aux[C].Nc)-1.5*log(fs[C].T)));
108 PetscScalar Evc = e*(Eqvc/e - kb*fs[C].T/e*(log(aux[C].Nv)-1.5*log(fs[C].T)));
109 PetscScalar Rc = mt->band->Recomb(pc,nc,fs[C].T);
110 PetscScalar nic = mt->band->nie(T);
111 PetscScalar phinc = Vc - log(fabs(nc)/nic)*Vt;
112 PetscScalar phipc = Vc + log(fabs(pc)/nic)*Vt;
113 PetscScalar qcc = -(-e*phinc-Eqcc)/(kb*T);
114 PetscScalar qvc = (-e*phipc-Eqvc)/(kb*T);
115
116 PetscScalar Ex = -((yb-yc)*Va + (yc-ya)*Vb +(ya-yb)*Vc)/(2*tri_area);
117 PetscScalar Ey = -((xc-xb)*Va + (xa-xc)*Vb +(xb-xa)*Vc)/(2*tri_area);
118 PetscScalar E = sqrt(Ex*Ex+Ey*Ey+1e-20);
119
120 PetscScalar G = 0;
121 PetscScalar IIn=0,IIp=0;
122 PetscScalar mun,mup;
123 PetscScalar Epnc=0,Etnc=0;
124 PetscScalar Eppc=0,Etpc=0;
125 PetscScalar Epna=0,Etna=0;
126 PetscScalar Eppa=0,Etpa=0;
127 PetscScalar Epnb=0,Etnb=0;
128 PetscScalar Eppb=0,Etpb=0;
129 PetscScalar Jna_norm,Jpa_norm;
130 PetscScalar Jnb_norm,Jpb_norm;
131 PetscScalar Jnc_norm,Jpc_norm;
132
133 #ifdef _FLUX1_
134 PetscScalar Jnc = In(Vt,Eca/e,Ecb/e,na,nb,ptri->edge_len[2]);
135 PetscScalar Jpc = Ip(Vt,Eva/e,Evb/e,pa,pb,ptri->edge_len[2]);
136 PetscScalar Jna = In(Vt,Ecb/e,Ecc/e,nb,nc,ptri->edge_len[0]);
137 PetscScalar Jpa = Ip(Vt,Evb/e,Evc/e,pb,pc,ptri->edge_len[0]);
138 PetscScalar Jnb = In(Vt,Ecc/e,Eca/e,nc,na,ptri->edge_len[1]);
139 PetscScalar Jpb = Ip(Vt,Evc/e,Eva/e,pc,pa,ptri->edge_len[1]);
140 PetscScalar Jn_scale = dmax(dmax(fabs(Jna),fabs(Jnb)),dmax(fabs(Jnc),nia*nia));
141 PetscScalar Jp_scale = dmax(dmax(fabs(Jpa),fabs(Jpb)),dmax(fabs(Jpc),nia*nia));
142 //cout<<Jn_scale<<" "<<Jnc<<" "<<Jna<<" "<<Jnb<<endl;
143 //cout<<Jp_scale<<" "<<Jpc<<" "<<Jpa<<" "<<Jpb<<endl;
144 Jnc /= Jn_scale;
145 Jpc /= Jp_scale;
146 Jna /= Jn_scale;
147 Jpa /= Jp_scale;
148 Jnb /= Jn_scale;
149 Jpb /= Jp_scale;
150 #endif
151 #ifdef _FLUX2_
152 PetscScalar Jnc = In(Vt,(Ecb-Eca)/e,na,nb,ptri->edge_len[2]);
153 PetscScalar Jpc = Ip(Vt,(Evb-Eva)/e,pa,pb,ptri->edge_len[2]);
154 PetscScalar Jna = In(Vt,(Ecc-Ecb)/e,nb,nc,ptri->edge_len[0]);
155 PetscScalar Jpa = Ip(Vt,(Evc-Evb)/e,pb,pc,ptri->edge_len[0]);
156 PetscScalar Jnb = In(Vt,(Eca-Ecc)/e,nc,na,ptri->edge_len[1]);
157 PetscScalar Jpb = Ip(Vt,(Eva-Evc)/e,pc,pa,ptri->edge_len[1]);
158 PetscScalar Jn_scale = dmax(dmax(fabs(Jna),fabs(Jnb)),dmax(fabs(Jnc),nia*nia));
159 PetscScalar Jp_scale = dmax(dmax(fabs(Jpa),fabs(Jpb)),dmax(fabs(Jpc),nia*nia));
160 Jnc /= Jn_scale;
161 Jpc /= Jp_scale;
162 Jna /= Jn_scale;
163 Jpa /= Jp_scale;
164 Jnb /= Jn_scale;
165 Jpb /= Jp_scale;
166 #endif
167
168 //flux along A-B
169 if(HighFieldMobility)
170 {
171 if(EJModel || IIType==EdotJ)
172 {
173 PetscScalar Jncx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jnc
174 + (Wca*Sax - Wca*Mb*Scx)*Jna
175 + (Wcb*Sbx - Wcb*Ma*Scx)*Jnb;
176 PetscScalar Jncy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jnc
177 + (Wca*Say - Wca*Mb*Scy)*Jna
178 + (Wcb*Sby - Wcb*Ma*Scy)*Jnb;
179 PetscScalar Jpcx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jpc
180 + (Wca*Sax - Wca*Mb*Scx)*Jpa
181 + (Wcb*Sbx - Wcb*Ma*Scx)*Jpb;
182 PetscScalar Jpcy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jpc
183 + (Wca*Say - Wca*Mb*Scy)*Jpa
184 + (Wcb*Sby - Wcb*Ma*Scy)*Jpb;
185 Jnc_norm = sqrt(Jncx*Jncx + Jncy*Jncy + 1e-100);
186 Jpc_norm = sqrt(Jpcx*Jpcx + Jpcy*Jpcy + 1e-100);
187 Epnc = (Ex*Jncx + Ey*Jncy)/Jnc_norm;
188 Etnc = (Ex*Jncy - Ey*Jncx)/Jnc_norm;
189 Eppc = (Ex*Jpcx + Ey*Jpcy)/Jpc_norm;
190 Etpc = (Ex*Jpcy - Ey*Jpcx)/Jpc_norm;
191 }
192 else
193 {
194 Epnc = fabs(Eca-Ecb)/e/ptri->edge_len[2]; //parallel electrical field for electron
195 Eppc = fabs(Eva-Evb)/e/ptri->edge_len[2]; //parallel electrical field for hole
196 //transvers electrical field for electron and hole
197 Etnc = fabs(Eca + ptri->edge_len[1]*cos(ptri->angle[0])*(Ecb-Eca)/ptri->edge_len[2] - Ecc)
198 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
199 Etpc = fabs(Eva + ptri->edge_len[1]*cos(ptri->angle[0])*(Evb-Eva)/ptri->edge_len[2] - Evc)
200 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
201 }
202
203 mt->mapping(&pzone->danode[A],&aux[A],0);
204 aux[A].mun = mt->mob->ElecMob(pa,na,T,dmax(0,Epnc),fabs(Etnc),T);
205 aux[A].mup = mt->mob->HoleMob(pa,na,T,dmax(0,Eppc),fabs(Etpc),T);
206
207 mt->mapping(&pzone->danode[B],&aux[B],0);
208 aux[B].mun = mt->mob->ElecMob(pb,nb,T,dmax(0,Epnc),fabs(Etnc),T);
209 aux[B].mup = mt->mob->HoleMob(pb,nb,T,dmax(0,Eppc),fabs(Etpc),T);
210 }
211
212 mun = 0.5*(aux[A].mun+aux[B].mun);
213 mup = 0.5*(aux[A].mup+aux[B].mup);
214
215
216 f[zofs[zone_index]+5*A+0] += 0.5*(aux[A].eps+aux[B].eps)*(Vb-Va)/ptri->edge_len[2]*ptri->d[2];
217 f[zofs[zone_index]+5*A+1] += mun*Jnc*Jn_scale*ptri->d[2];
218 f[zofs[zone_index]+5*A+2] += - mup*Jpc*Jp_scale*ptri->d[2];
219 f[zofs[zone_index]+5*A+3] += - QNFactor*ptri->d[2]*bn*qV(qca,qcb,ptri->edge_len[2]);
220 f[zofs[zone_index]+5*A+4] += QPFactor*ptri->d[2]*bp*qV(qva,qvb,ptri->edge_len[2]);
221
222 f[zofs[zone_index]+5*B+0] += - 0.5*(aux[A].eps+aux[B].eps)*(Vb-Va)/ptri->edge_len[2]*ptri->d[2];
223 f[zofs[zone_index]+5*B+1] += - mun*Jnc*Jn_scale*ptri->d[2];
224 f[zofs[zone_index]+5*B+2] += mup*Jpc*Jp_scale*ptri->d[2];
225 f[zofs[zone_index]+5*B+3] += - QNFactor*ptri->d[2]*bn*qV(qcb,qca,ptri->edge_len[2]);
226 f[zofs[zone_index]+5*B+4] += QPFactor*ptri->d[2]*bp*qV(qvb,qva,ptri->edge_len[2]);
227
228 //flux along B-C
229 if(HighFieldMobility)
230 {
231 if(EJModel || IIType==EdotJ)
232 {
233 PetscScalar Jnax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jna
234 + (Wab*Sbx - Wab*Mc*Sax)*Jnb
235 + (Wac*Scx - Wac*Mb*Sax)*Jnc;
236 PetscScalar Jnay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jna
237 + (Wab*Sby - Wab*Mc*Say)*Jnb
238 + (Wac*Scy - Wac*Mb*Say)*Jnc;
239 PetscScalar Jpax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jpa
240 + (Wab*Sbx - Wab*Mc*Sax)*Jpb
241 + (Wac*Scx - Wac*Mb*Sax)*Jpc;
242 PetscScalar Jpay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jpa
243 + (Wab*Sby - Wab*Mc*Say)*Jpb
244 + (Wac*Scy - Wac*Mb*Say)*Jpc;
245 Jna_norm = sqrt(Jnax*Jnax + Jnay*Jnay + 1e-100);
246 Jpa_norm = sqrt(Jpax*Jpax + Jpay*Jpay + 1e-100);
247 Epna = (Ex*Jnax + Ey*Jnay)/Jna_norm;
248 Etna = (Ex*Jnay - Ey*Jnax)/Jna_norm;
249 Eppa = (Ex*Jpax + Ey*Jpay)/Jpa_norm;
250 Etpa = (Ex*Jpay - Ey*Jpax)/Jpa_norm;
251 }
252 else
253 {
254 Epna = fabs(Ecb-Ecc)/e/ptri->edge_len[0]; //parallel electrical field for electron
255 Eppa = fabs(Evb-Evc)/e/ptri->edge_len[0]; //parallel electrical field for hole
256 //transvers electrical field for electron and hole
257 Etna = fabs(Ecb + ptri->edge_len[2]*cos(ptri->angle[1])*(Ecc-Ecb)/ptri->edge_len[0] - Eca)
258 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
259 Etpa = fabs(Evb + ptri->edge_len[2]*cos(ptri->angle[1])*(Evc-Evb)/ptri->edge_len[0] - Eva)
260 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
261 }
262 mt->mapping(&pzone->danode[B],&aux[B],0);
263 aux[B].mun = mt->mob->ElecMob(pb,nb,T,dmax(0,Epna),fabs(Etna),T);
264 aux[B].mup = mt->mob->HoleMob(pb,nb,T,dmax(0,Eppa),fabs(Etpa),T);
265
266 mt->mapping(&pzone->danode[C],&aux[C],0);
267 aux[C].mun = mt->mob->ElecMob(pc,nc,T,dmax(0,Epna),fabs(Etna),T);
268 aux[C].mup = mt->mob->HoleMob(pc,nc,T,dmax(0,Eppa),fabs(Etpa),T);
269 }
270 mun = 0.5*(aux[B].mun+aux[C].mun);
271 mup = 0.5*(aux[B].mup+aux[C].mup);
272
273
274 f[zofs[zone_index]+5*B+0] += 0.5*(aux[B].eps+aux[C].eps)*(Vc-Vb)/ptri->edge_len[0]*ptri->d[0];
275 f[zofs[zone_index]+5*B+1] += mun*Jna*Jn_scale*ptri->d[0];
276 f[zofs[zone_index]+5*B+2] += - mup*Jpa*Jp_scale*ptri->d[0];
277 f[zofs[zone_index]+5*B+3] += - QNFactor*ptri->d[0]*bn*qV(qcb,qcc,ptri->edge_len[0]);
278 f[zofs[zone_index]+5*B+4] += QPFactor*ptri->d[0]*bp*qV(qvb,qvc,ptri->edge_len[0]);
279
280 f[zofs[zone_index]+5*C+0] += - 0.5*(aux[B].eps+aux[C].eps)*(Vc-Vb)/ptri->edge_len[0]*ptri->d[0];
281 f[zofs[zone_index]+5*C+1] += - mun*Jna*Jn_scale*ptri->d[0];
282 f[zofs[zone_index]+5*C+2] += mup*Jpa*Jp_scale*ptri->d[0];
283 f[zofs[zone_index]+5*C+3] += - QNFactor*ptri->d[0]*bn*qV(qcc,qcb,ptri->edge_len[0]);
284 f[zofs[zone_index]+5*C+4] += QPFactor*ptri->d[0]*bp*qV(qvc,qvb,ptri->edge_len[0]);
285
286 //flux along C-A
287 if(HighFieldMobility)
288 {
289 if(EJModel || IIType==EdotJ)
290 {
291 PetscScalar Jnbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jnb
292 + (Wbc*Scx - Wbc*Ma*Sbx)*Jnc
293 + (Wba*Sax - Wba*Mc*Sbx)*Jna;
294 PetscScalar Jnby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jnb
295 + (Wbc*Scy - Wbc*Ma*Sby)*Jnc
296 + (Wba*Say - Wba*Mc*Sby)*Jna;
297 PetscScalar Jpbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jpb
298 + (Wbc*Scx - Wbc*Ma*Sbx)*Jpc
299 + (Wba*Sax - Wba*Mc*Sbx)*Jpa;
300 PetscScalar Jpby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jpb
301 + (Wbc*Scy - Wbc*Ma*Sby)*Jpc
302 + (Wba*Say - Wba*Mc*Sby)*Jpa;
303 Jnb_norm = sqrt(Jnbx*Jnbx + Jnby*Jnby + 1e-100);
304 Jpb_norm = sqrt(Jpbx*Jpbx + Jpby*Jpby + 1e-100);
305 Epnb = (Ex*Jnbx + Ey*Jnby)/Jnb_norm;
306 Etnb = (Ex*Jnby - Ey*Jnbx)/Jnb_norm;
307 Eppb = (Ex*Jpbx + Ey*Jpby)/Jpb_norm;
308 Etpb = (Ex*Jpby - Ey*Jpbx)/Jpb_norm;
309
310 }
311 else
312 {
313 Epnb = fabs(Ecc-Eca)/e/ptri->edge_len[1]; //parallel electrical field for electron
314 Eppb = fabs(Evc-Eva)/e/ptri->edge_len[1]; //parallel electrical field for hole
315 //transvers electrical field for electron and hole
316 Etnb = fabs(Ecc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eca-Ecc)/ptri->edge_len[1] - Ecb)
317 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
318 Etpb = fabs(Evc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eva-Evc)/ptri->edge_len[1] - Evb)
319 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
320 }
321 mt->mapping(&pzone->danode[C],&aux[C],0);
322 aux[C].mun = mt->mob->ElecMob(pc,nc,T,dmax(0,Epnb),fabs(Etnb),T);
323 aux[C].mup = mt->mob->HoleMob(pc,nc,T,dmax(0,Eppb),fabs(Etpb),T);
324
325 mt->mapping(&pzone->danode[A],&aux[A],0);
326 aux[A].mun = mt->mob->ElecMob(pa,na,T,dmax(0,Epnb),fabs(Etnb),T);
327 aux[A].mup = mt->mob->HoleMob(pa,na,T,dmax(0,Eppb),fabs(Etpb),T);
328 }
329 mun = 0.5*(aux[C].mun+aux[A].mun);
330 mup = 0.5*(aux[C].mup+aux[A].mup);
331
332
333 f[zofs[zone_index]+5*C+0] += 0.5*(aux[C].eps+aux[A].eps)*(Va-Vc)/ptri->edge_len[1]*ptri->d[1];
334 f[zofs[zone_index]+5*C+1] += mun*Jnb*Jn_scale*ptri->d[1];
335 f[zofs[zone_index]+5*C+2] += - mup*Jpb*Jp_scale*ptri->d[1];
336 f[zofs[zone_index]+5*C+3] += - QNFactor*ptri->d[1]*bn*qV(qcc,qca,ptri->edge_len[1]);
337 f[zofs[zone_index]+5*C+4] += QPFactor*ptri->d[1]*bp*qV(qvc,qva,ptri->edge_len[1]);
338
339 f[zofs[zone_index]+5*A+0] += - 0.5*(aux[C].eps+aux[A].eps)*(Va-Vc)/ptri->edge_len[1]*ptri->d[1];
340 f[zofs[zone_index]+5*A+1] += - mun*Jnb*Jn_scale*ptri->d[1];
341 f[zofs[zone_index]+5*A+2] += mup*Jpb*Jp_scale*ptri->d[1];
342 f[zofs[zone_index]+5*A+3] += - QNFactor*ptri->d[1]*bn*qV(qca,qcc,ptri->edge_len[1]);
343 f[zofs[zone_index]+5*A+4] += QPFactor*ptri->d[1]*bp*qV(qva,qvc,ptri->edge_len[1]);
344
345 //---------------------------------------------------------------------------
346 // Recombination item
347 //---------------------------------------------------------------------------
348
349 S = 0.25*ptri->d[2]*ptri->edge_len[2] + 0.25*ptri->d[1]*ptri->edge_len[1];
350 f[zofs[zone_index]+5*A+1] += -Ra*S;
351 f[zofs[zone_index]+5*A+2] += -Ra*S;
352
353 S = 0.25*ptri->d[0]*ptri->edge_len[0] + 0.25*ptri->d[2]*ptri->edge_len[2];
354 f[zofs[zone_index]+5*B+1] += -Rb*S;
355 f[zofs[zone_index]+5*B+2] += -Rb*S;
356
357 S = 0.25*ptri->d[1]*ptri->edge_len[1] + 0.25*ptri->d[0]*ptri->edge_len[0];
358 f[zofs[zone_index]+5*C+1] += -Rc*S;
359 f[zofs[zone_index]+5*C+2] += -Rc*S;
360 }
361
362
363
364 void SMCZone::J1Q_Tri_qddm(Tri *ptri,PetscScalar *x,Mat *jtmp, vector<int> & zofs)
365 {
366 int z = zone_index;
367 int A = ptri->node[0];
368 int B = ptri->node[1];
369 int C = ptri->node[2];
370
371 PetscScalar xa = pzone->danode[ptri->node[0]].x;
372 PetscScalar xb = pzone->danode[ptri->node[1]].x;
373 PetscScalar xc = pzone->danode[ptri->node[2]].x;
374 PetscScalar ya = pzone->danode[ptri->node[0]].y;
375 PetscScalar yb = pzone->danode[ptri->node[1]].y;
376 PetscScalar yc = pzone->danode[ptri->node[2]].y;
377 PetscScalar tri_area = ptri->area;
378 PetscScalar Sax = (xc-xb)/ptri->edge_len[0];
379 PetscScalar Say = (yc-yb)/ptri->edge_len[0];
380 PetscScalar Sbx = (xa-xc)/ptri->edge_len[1];
381 PetscScalar Sby = (ya-yc)/ptri->edge_len[1];
382 PetscScalar Scx = (xb-xa)/ptri->edge_len[2];
383 PetscScalar Scy = (yb-ya)/ptri->edge_len[2];
384 PetscScalar Wab = ptri->d[1]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
385 PetscScalar Wac = ptri->d[2]/((ptri->d[1]+ptri->d[2])*sin(ptri->angle[1])*sin(ptri->angle[1]));
386 PetscScalar Wbc = ptri->d[2]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[0])*sin(ptri->angle[0]));
387 PetscScalar Wba = ptri->d[0]/((ptri->d[0]+ptri->d[2])*sin(ptri->angle[2])*sin(ptri->angle[2]));
388 PetscScalar Wca = ptri->d[0]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[1])*sin(ptri->angle[1]));
389 PetscScalar Wcb = ptri->d[1]/((ptri->d[0]+ptri->d[1])*sin(ptri->angle[0])*sin(ptri->angle[0]));
390 PetscScalar Ma = -cos(ptri->angle[0]);
391 PetscScalar Mb = -cos(ptri->angle[1]);
392 PetscScalar Mc = -cos(ptri->angle[2]);
393
394 Mat5 J1,J2,J3;
395 Vec5 scale;
396 PetscInt index1[5],index2[5],index3[5];
397 index1[0] = zofs[z]+5*A+0;
398 index1[1] = zofs[z]+5*A+1;
399 index1[2] = zofs[z]+5*A+2;
400 index1[3] = zofs[z]+5*A+3;
401 index1[4] = zofs[z]+5*A+4;
402
403 index2[0] = zofs[z]+5*B+0;
404 index2[1] = zofs[z]+5*B+1;
405 index2[2] = zofs[z]+5*B+2;
406 index2[3] = zofs[z]+5*B+3;
407 index2[4] = zofs[z]+5*B+4;
408
409 index3[0] = zofs[z]+5*C+0;
410 index3[1] = zofs[z]+5*C+1;
411 index3[2] = zofs[z]+5*C+2;
412 index3[3] = zofs[z]+5*C+3;
413 index3[4] = zofs[z]+5*C+4;
414
415 //---------------------------------------------------------------------------
416
417 PetscScalar kb = mt->kb;
418 PetscScalar hbar = mt->hbar;
419 PetscScalar e = mt->e;
420 PetscScalar T = (fs[A].T+fs[B].T+fs[C].T)/3.0;
421 PetscScalar Eg = mt->band->Eg(T);
422 PetscScalar Vt = kb*T/e;
423 PetscScalar me = mt->band->EffecElecMass(T);
424 PetscScalar mh = mt->band->EffecHoleMass(T);
425 PetscScalar gamman = mt->band->Gamman();
426 PetscScalar gammap = mt->band->Gammap();
427 PetscScalar bn = gamman*hbar*hbar/(6*e*me);
428 PetscScalar bp = gammap*hbar*hbar/(6*e*mh);
429 PetscScalar S;
430
431 //the indepedent variable number
432 adtl::AutoDScalar::numdir=15;
433 //synchronize with material database
434 mt->set_ad_num(adtl::AutoDScalar::numdir);
435 AutoDScalar TD = T; // dummy ad variable.
436 AutoDScalar F[5];
437
438 AutoDScalar Va = x[zofs[zone_index]+5*A+0]; //potential of node A
439 AutoDScalar na = x[zofs[zone_index]+5*A+1]; //electron density of node A
440 AutoDScalar pa = x[zofs[zone_index]+5*A+2]; //hole density of node A
441 AutoDScalar Eqca = x[zofs[zone_index]+5*A+3]; //quantum conduction band of node A
442 AutoDScalar Eqva = x[zofs[zone_index]+5*A+4]; //quantum valence band of node A
443 Va.setADValue(0,1.0);
444 na.setADValue(1,1.0);
445 pa.setADValue(2,1.0);
446 Eqca.setADValue(3,1.0);
447 Eqva.setADValue(4,1.0);
448 mt->mapping(&pzone->danode[A],&aux[A],0);
449 AutoDScalar Eca = e*(Eqca/e + kb*fs[A].T/e*(log(aux[A].Nc)-1.5*log(fs[A].T)));
450 AutoDScalar Eva = e*(Eqva/e - kb*fs[A].T/e*(log(aux[A].Nv)-1.5*log(fs[A].T)));
451 AutoDScalar Ra = mt->band->Recomb(pa,na,fs[A].T);
452 PetscScalar nia = mt->band->nie(T);
453 AutoDScalar phina = Va - log(fabs(na)/nia)*Vt;
454 AutoDScalar phipa = Va + log(fabs(pa)/nia)*Vt;
455 AutoDScalar qca = -(-e*phina-Eqca)/(kb*T);
456 AutoDScalar qva = (-e*phipa-Eqva)/(kb*T);
457
458 AutoDScalar Vb = x[zofs[zone_index]+5*B+0]; //potential of node B
459 AutoDScalar nb = x[zofs[zone_index]+5*B+1]; //electron density of node B
460 AutoDScalar pb = x[zofs[zone_index]+5*B+2]; //hole density of node B
461 AutoDScalar Eqcb = x[zofs[zone_index]+5*B+3]; //quantum conduction band of node B
462 AutoDScalar Eqvb = x[zofs[zone_index]+5*B+4]; //quantum valence band of node B
463 Vb.setADValue(5,1.0);
464 nb.setADValue(6,1.0);
465 pb.setADValue(7,1.0);
466 Eqcb.setADValue(8,1.0);
467 Eqvb.setADValue(9,1.0);
468 mt->mapping(&pzone->danode[B],&aux[B],0);
469 AutoDScalar Ecb = e*(Eqcb/e + kb*fs[B].T/e*(log(aux[B].Nc)-1.5*log(fs[B].T)));
470 AutoDScalar Evb = e*(Eqvb/e - kb*fs[B].T/e*(log(aux[B].Nv)-1.5*log(fs[B].T)));
471 AutoDScalar Rb = mt->band->Recomb(pb,nb,fs[B].T);
472 PetscScalar nib = mt->band->nie(T);
473 AutoDScalar phinb = Vb - log(fabs(nb)/nib)*Vt;
474 AutoDScalar phipb = Vb + log(fabs(pb)/nib)*Vt;
475 AutoDScalar qcb = -(-e*phinb-Eqcb)/(kb*T);
476 AutoDScalar qvb = (-e*phipb-Eqvb)/(kb*T);
477
478 AutoDScalar Vc = x[zofs[zone_index]+5*C+0]; //potential of node C
479 AutoDScalar nc = x[zofs[zone_index]+5*C+1]; //electron density of node C
480 AutoDScalar pc = x[zofs[zone_index]+5*C+2]; //hole density of node C
481 AutoDScalar Eqcc = x[zofs[zone_index]+5*C+3]; //quantum conduction band of node C
482 AutoDScalar Eqvc = x[zofs[zone_index]+5*C+4]; //quantum valence band of node C
483 Vc.setADValue(10,1.0);
484 nc.setADValue(11,1.0);
485 pc.setADValue(12,1.0);
486 Eqcc.setADValue(13,1.0);
487 Eqvc.setADValue(14,1.0);
488 mt->mapping(&pzone->danode[C],&aux[C],0);
489 AutoDScalar Ecc = e*(Eqcc/e + kb*fs[C].T/e*(log(aux[C].Nc)-1.5*log(fs[C].T)));
490 AutoDScalar Evc = e*(Eqvc/e - kb*fs[C].T/e*(log(aux[C].Nv)-1.5*log(fs[C].T)));
491 AutoDScalar Rc = mt->band->Recomb(pc,nc,fs[C].T);
492 PetscScalar nic = mt->band->nie(T);
493 AutoDScalar phinc = Vc - log(fabs(nc)/nic)*Vt;
494 AutoDScalar phipc = Vc + log(fabs(pc)/nic)*Vt;
495 AutoDScalar qcc = -(-e*phinc-Eqcc)/(kb*T);
496 AutoDScalar qvc = (-e*phipc-Eqvc)/(kb*T);
497
498 AutoDScalar Ex = -((yb-yc)*Va + (yc-ya)*Vb +(ya-yb)*Vc)/(2*tri_area);
499 AutoDScalar Ey = -((xc-xb)*Va + (xa-xc)*Vb +(xb-xa)*Vc)/(2*tri_area);
500 AutoDScalar E = sqrt(Ex*Ex+Ey*Ey+1e-20);
501
502 AutoDScalar G = 0;
503 AutoDScalar IIn=0,IIp=0;
504 AutoDScalar mun,mup;
505 AutoDScalar Epnc=0,Etnc=0;
506 AutoDScalar Eppc=0,Etpc=0;
507 AutoDScalar Epna=0,Etna=0;
508 AutoDScalar Eppa=0,Etpa=0;
509 AutoDScalar Epnb=0,Etnb=0;
510 AutoDScalar Eppb=0,Etpb=0;
511 AutoDScalar Jna_norm=0,Jpa_norm=0;
512 AutoDScalar Jnb_norm=0,Jpb_norm=0;
513 AutoDScalar Jnc_norm=0,Jpc_norm=0;
514
515 #ifdef _FLUX1_
516 AutoDScalar Jnc = In(Vt,Eca/e,Ecb/e,na,nb,ptri->edge_len[2]);
517 AutoDScalar Jpc = Ip(Vt,Eva/e,Evb/e,pa,pb,ptri->edge_len[2]);
518 AutoDScalar Jna = In(Vt,Ecb/e,Ecc/e,nb,nc,ptri->edge_len[0]);
519 AutoDScalar Jpa = Ip(Vt,Evb/e,Evc/e,pb,pc,ptri->edge_len[0]);
520 AutoDScalar Jnb = In(Vt,Ecc/e,Eca/e,nc,na,ptri->edge_len[1]);
521 AutoDScalar Jpb = Ip(Vt,Evc/e,Eva/e,pc,pa,ptri->edge_len[1]);
522 PetscScalar Jn_scale = adtl::fmax(adtl::fmax(fabs(Jna.getValue()),fabs(Jnb.getValue())),
523 adtl::fmax(fabs(Jnc.getValue()),nia*nia));
524 PetscScalar Jp_scale = adtl::fmax(adtl::fmax(fabs(Jpa.getValue()),fabs(Jpb.getValue())),
525 adtl::fmax(fabs(Jpc.getValue()),nia*nia));
526
527 Jnc = Jnc/Jn_scale;
528 Jpc = Jpc/Jp_scale;
529 Jna = Jna/Jn_scale;
530 Jpa = Jpa/Jp_scale;
531 Jnb = Jnb/Jn_scale;
532 Jpb = Jpb/Jp_scale;
533
534 #endif
535 #ifdef _FLUX2_
536 AutoDScalar Jnc = In(Vt,(Ecb-Eca)/e,na,nb,ptri->edge_len[2]);
537 AutoDScalar Jpc = Ip(Vt,(Evb-Eva)/e,pa,pb,ptri->edge_len[2]);
538 AutoDScalar Jna = In(Vt,(Ecc-Ecb)/e,nb,nc,ptri->edge_len[0]);
539 AutoDScalar Jpa = Ip(Vt,(Evc-Evb)/e,pb,pc,ptri->edge_len[0]);
540 AutoDScalar Jnb = In(Vt,(Eca-Ecc)/e,nc,na,ptri->edge_len[1]);
541 AutoDScalar Jpb = Ip(Vt,(Eva-Evc)/e,pc,pa,ptri->edge_len[1]);
542 PetscScalar Jn_scale = adtl::fmax(adtl::fmax(fabs(Jna.getValue()),fabs(Jnb.getValue())),
543 adtl::fmax(fabs(Jnc.getValue()),nia*nia));
544 PetscScalar Jp_scale = adtl::fmax(adtl::fmax(fabs(Jpa.getValue()),fabs(Jpb.getValue())),
545 adtl::fmax(fabs(Jpc.getValue()),nia*nia));
546 Jnc /= Jn_scale;
547 Jpc /= Jp_scale;
548 Jna /= Jn_scale;
549 Jpa /= Jp_scale;
550 Jnb /= Jn_scale;
551 Jpb /= Jp_scale;
552 #endif
553
554 //---------------------------------------------------------------------------
555 //flux along A-B
556 //---------------------------------------------------------------------------
557 if(HighFieldMobility)
558 {
559 if(EJModel || IIType==EdotJ)
560 {
561 AutoDScalar Jncx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jnc
562 + (Wca*Sax - Wca*Mb*Scx)*Jna
563 + (Wcb*Sbx - Wcb*Ma*Scx)*Jnb;
564 AutoDScalar Jncy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jnc
565 + (Wca*Say - Wca*Mb*Scy)*Jna
566 + (Wcb*Sby - Wcb*Ma*Scy)*Jnb;
567 AutoDScalar Jpcx = ((Wca+Wcb)*Scx - Wca*Mb*Sax - Wcb*Ma*Sbx)*Jpc
568 + (Wca*Sax - Wca*Mb*Scx)*Jpa
569 + (Wcb*Sbx - Wcb*Ma*Scx)*Jpb;
570 AutoDScalar Jpcy = ((Wca+Wcb)*Scy - Wca*Mb*Say - Wcb*Ma*Sby)*Jpc
571 + (Wca*Say - Wca*Mb*Scy)*Jpa
572 + (Wcb*Sby - Wcb*Ma*Scy)*Jpb;
573 Jnc_norm = sqrt(Jncx*Jncx + Jncy*Jncy + 1e-100);
574 Jpc_norm = sqrt(Jpcx*Jpcx + Jpcy*Jpcy + 1e-100);
575 Epnc = Ex*(Jncx/Jnc_norm) + Ey*(Jncy/Jnc_norm);
576 Etnc = Ex*(Jncy/Jnc_norm) - Ey*(Jncx/Jnc_norm);
577 Eppc = Ex*(Jpcx/Jpc_norm) + Ey*(Jpcy/Jpc_norm);
578 Etpc = Ex*(Jpcy/Jpc_norm) - Ey*(Jpcx/Jpc_norm);
579 }
580 else
581 {
582 Epnc = fabs(Eca-Ecb)/e/ptri->edge_len[2]; //parallel electrical field for electron
583 Eppc = fabs(Eva-Evb)/e/ptri->edge_len[2]; //parallel electrical field for hole
584 //transvers electrical field for electron and hole
585 Etnc = fabs(Eca + ptri->edge_len[1]*cos(ptri->angle[0])*(Ecb-Eca)/ptri->edge_len[2] - Ecc)
586 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
587 Etpc = fabs(Eva + ptri->edge_len[1]*cos(ptri->angle[0])*(Evb-Eva)/ptri->edge_len[2] - Evc)
588 /(ptri->edge_len[1]*sin(ptri->angle[0]))/e;
589 }
590 mt->mapping(&pzone->danode[A],&aux[A],0);
591 AutoDScalar munCA = mt->mob->ElecMob(pa,na,TD,adtl::fmax(0,Epnc),fabs(Etnc),TD);
592 AutoDScalar mupCA = mt->mob->HoleMob(pa,na,TD,adtl::fmax(0,Eppc),fabs(Etpc),TD);
593
594 mt->mapping(&pzone->danode[B],&aux[B],0);
595 AutoDScalar munCB = mt->mob->ElecMob(pb,nb,TD,adtl::fmax(0,Epnc),fabs(Etnc),TD);
596 AutoDScalar mupCB = mt->mob->HoleMob(pb,nb,TD,adtl::fmax(0,Eppc),fabs(Etpc),TD);
597
598 mun = 0.5*(munCA+munCB);
599 mup = 0.5*(mupCA+mupCB);
600 }
601 else
602 {
603 mun = 0.5*(aux[A].mun+aux[B].mun);
604 mup = 0.5*(aux[A].mup+aux[B].mup);
605 }
606
607 //---------------------------------------------------------------------------
608 F[0] = 0.5*(aux[A].eps+aux[B].eps)*(Vb-Va)/ptri->edge_len[2]*ptri->d[2];
609 F[1] = mun*Jnc*Jn_scale*ptri->d[2];
610 F[2] = - mup*Jpc*Jp_scale*ptri->d[2];
611 F[3] = - QNFactor*ptri->d[2]*bn*qV(qca,qcb,ptri->edge_len[2]);
612 F[4] = QPFactor*ptri->d[2]*bp*qV(qva,qvb,ptri->edge_len[2]);
613
614 for(int i=0;i<5;i++)
615 for(int j=0;j<5;j++)
616 {
617 J1.m[5*i+j] = F[i].getADValue(0+j);
618 J2.m[5*i+j] = F[i].getADValue(5+j);
619 J3.m[5*i+j] = F[i].getADValue(10+j);
620 }
621 MatSetValues(*jtmp,5,index1,5,index1,J1.m,ADD_VALUES);
622 MatSetValues(*jtmp,5,index1,5,index2,J2.m,ADD_VALUES);
623 MatSetValues(*jtmp,5,index1,5,index3,J3.m,ADD_VALUES);
624
625 F[0] = - 0.5*(aux[A].eps+aux[B].eps)*(Vb-Va)/ptri->edge_len[2]*ptri->d[2];
626 F[1] = - mun*Jnc*Jn_scale*ptri->d[2];
627 F[2] = mup*Jpc*Jp_scale*ptri->d[2];
628 F[3] = - QNFactor*ptri->d[2]*bn*qV(qcb,qca,ptri->edge_len[2]);
629 F[4] = QPFactor*ptri->d[2]*bp*qV(qvb,qva,ptri->edge_len[2]);
630 for(int i=0;i<5;i++)
631 for(int j=0;j<5;j++)
632 {
633 J1.m[5*i+j] = F[i].getADValue(0+j);
634 J2.m[5*i+j] = F[i].getADValue(5+j);
635 J3.m[5*i+j] = F[i].getADValue(10+j);
636 }
637 MatSetValues(*jtmp,5,index2,5,index1,J1.m,ADD_VALUES);
638 MatSetValues(*jtmp,5,index2,5,index2,J2.m,ADD_VALUES);
639 MatSetValues(*jtmp,5,index2,5,index3,J3.m,ADD_VALUES);
640
641 //---------------------------------------------------------------------------
642 //flux along B-C
643 //---------------------------------------------------------------------------
644 if(HighFieldMobility)
645 {
646 if(EJModel || IIType==EdotJ)
647 {
648 AutoDScalar Jnax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jna
649 + (Wab*Sbx - Wab*Mc*Sax)*Jnb
650 + (Wac*Scx - Wac*Mb*Sax)*Jnc;
651 AutoDScalar Jnay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jna
652 + (Wab*Sby - Wab*Mc*Say)*Jnb
653 + (Wac*Scy - Wac*Mb*Say)*Jnc;
654 AutoDScalar Jpax = ((Wab+Wac)*Sax - Wab*Mc*Sbx - Wac*Mb*Scx)*Jpa
655 + (Wab*Sbx - Wab*Mc*Sax)*Jpb
656 + (Wac*Scx - Wac*Mb*Sax)*Jpc;
657 AutoDScalar Jpay = ((Wab+Wac)*Say - Wab*Mc*Sby - Wac*Mb*Scy)*Jpa
658 + (Wab*Sby - Wab*Mc*Say)*Jpb
659 + (Wac*Scy - Wac*Mb*Say)*Jpc;
660 Jna_norm = sqrt(Jnax*Jnax + Jnay*Jnay + 1e-100);
661 Jpa_norm = sqrt(Jpax*Jpax + Jpay*Jpay + 1e-100);
662 Epna = Ex*(Jnax/Jna_norm) + Ey*(Jnay/Jna_norm);
663 Etna = Ex*(Jnay/Jna_norm) - Ey*(Jnax/Jna_norm);
664 Eppa = Ex*(Jpax/Jpa_norm) + Ey*(Jpay/Jpa_norm);
665 Etpa = Ex*(Jpay/Jpa_norm) - Ey*(Jpax/Jpa_norm);
666 }
667 else
668 {
669 Epna = fabs(Ecb-Ecc)/e/ptri->edge_len[0]; //parallel electrical field for electron
670 Eppa = fabs(Evb-Evc)/e/ptri->edge_len[0]; //parallel electrical field for hole
671 //transvers electrical field for electron and hole
672 Etna = fabs(Ecb + ptri->edge_len[2]*cos(ptri->angle[1])*(Ecc-Ecb)/ptri->edge_len[0] - Eca)
673 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
674 Etpa = fabs(Evb + ptri->edge_len[2]*cos(ptri->angle[1])*(Evc-Evb)/ptri->edge_len[0] - Eva)
675 /(ptri->edge_len[2]*sin(ptri->angle[1]))/e;
676 }
677 mt->mapping(&pzone->danode[B],&aux[B],0);
678 AutoDScalar munAB = mt->mob->ElecMob(pb,nb,TD,adtl::fmax(0,Epna),fabs(Etna),TD);
679 AutoDScalar mupAB = mt->mob->HoleMob(pb,nb,TD,adtl::fmax(0,Eppa),fabs(Etpa),TD);
680
681 mt->mapping(&pzone->danode[C],&aux[C],0);
682 AutoDScalar munAC = mt->mob->ElecMob(pc,nc,TD,adtl::fmax(0,Epna),fabs(Etna),TD);
683 AutoDScalar mupAC = mt->mob->HoleMob(pc,nc,TD,adtl::fmax(0,Eppa),fabs(Etpa),TD);
684
685 mun = 0.5*(munAB+munAC);
686 mup = 0.5*(mupAB+mupAC);
687 }
688 else
689 {
690 mun = 0.5*(aux[B].mun+aux[C].mun);
691 mup = 0.5*(aux[B].mup+aux[C].mup);
692 }
693 F[0] = 0.5*(aux[B].eps+aux[C].eps)*(Vc-Vb)/ptri->edge_len[0]*ptri->d[0];
694 F[1] = mun*Jna*Jn_scale*ptri->d[0];
695 F[2] = - mup*Jpa*Jp_scale*ptri->d[0];
696 F[3] = - QNFactor*ptri->d[0]*bn*qV(qcb,qcc,ptri->edge_len[0]);
697 F[4] = QPFactor*ptri->d[0]*bp*qV(qvb,qvc,ptri->edge_len[0]);
698 for(int i=0;i<5;i++)
699 for(int j=0;j<5;j++)
700 {
701 J1.m[5*i+j] = F[i].getADValue(0+j);
702 J2.m[5*i+j] = F[i].getADValue(5+j);
703 J3.m[5*i+j] = F[i].getADValue(10+j);
704 }
705 MatSetValues(*jtmp,5,index2,5,index1,J1.m,ADD_VALUES);
706 MatSetValues(*jtmp,5,index2,5,index2,J2.m,ADD_VALUES);
707 MatSetValues(*jtmp,5,index2,5,index3,J3.m,ADD_VALUES);
708
709 F[0] = - 0.5*(aux[B].eps+aux[C].eps)*(Vc-Vb)/ptri->edge_len[0]*ptri->d[0];
710 F[1] = - mun*Jna*Jn_scale*ptri->d[0];
711 F[2] = mup*Jpa*Jp_scale*ptri->d[0];
712 F[3] = - QNFactor*ptri->d[0]*bn*qV(qcc,qcb,ptri->edge_len[0]);
713 F[4] = QPFactor*ptri->d[0]*bp*qV(qvc,qvb,ptri->edge_len[0]);
714 for(int i=0;i<5;i++)
715 for(int j=0;j<5;j++)
716 {
717 J1.m[5*i+j] = F[i].getADValue(0+j);
718 J2.m[5*i+j] = F[i].getADValue(5+j);
719 J3.m[5*i+j] = F[i].getADValue(10+j);
720 }
721 MatSetValues(*jtmp,5,index3,5,index1,J1.m,ADD_VALUES);
722 MatSetValues(*jtmp,5,index3,5,index2,J2.m,ADD_VALUES);
723 MatSetValues(*jtmp,5,index3,5,index3,J3.m,ADD_VALUES);
724
725 //---------------------------------------------------------------------------
726 //flux along C-A
727 //---------------------------------------------------------------------------
728 if(HighFieldMobility)
729 {
730 if(EJModel || IIType==EdotJ)
731 {
732 AutoDScalar Jnbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jnb
733 + (Wbc*Scx - Wbc*Ma*Sbx)*Jnc
734 + (Wba*Sax - Wba*Mc*Sbx)*Jna;
735 AutoDScalar Jnby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jnb
736 + (Wbc*Scy - Wbc*Ma*Sby)*Jnc
737 + (Wba*Say - Wba*Mc*Sby)*Jna;
738 AutoDScalar Jpbx = ((Wbc+Wba)*Sbx - Wbc*Ma*Scx - Wba*Mc*Sax)*Jpb
739 + (Wbc*Scx - Wbc*Ma*Sbx)*Jpc
740 + (Wba*Sax - Wba*Mc*Sbx)*Jpa;
741 AutoDScalar Jpby = ((Wbc+Wba)*Sby - Wbc*Ma*Scy - Wba*Mc*Say)*Jpb
742 + (Wbc*Scy - Wbc*Ma*Sby)*Jpc
743 + (Wba*Say - Wba*Mc*Sby)*Jpa;
744 Jnb_norm = sqrt(Jnbx*Jnbx + Jnby*Jnby + 1e-100);
745 Jpb_norm = sqrt(Jpbx*Jpbx + Jpby*Jpby + 1e-100);
746 Epnb = Ex*(Jnbx/Jnb_norm) + Ey*(Jnby/Jnb_norm);
747 Etnb = Ex*(Jnby/Jnb_norm) - Ey*(Jnbx/Jnb_norm);
748 Eppb = Ex*(Jpbx/Jpb_norm) + Ey*(Jpby/Jpb_norm);
749 Etpb = Ex*(Jpby/Jpb_norm) - Ey*(Jpbx/Jpb_norm);
750 }
751 else
752 {
753 Epnb = fabs(Ecc-Eca)/e/ptri->edge_len[1]; //parallel electrical field for electron
754 Eppb = fabs(Evc-Eva)/e/ptri->edge_len[1]; //parallel electrical field for hole
755 //transvers electrical field for electron and hole
756 Etnb = fabs(Ecc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eca-Ecc)/ptri->edge_len[1] - Ecb)
757 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
758 Etpb = fabs(Evc + ptri->edge_len[0]*cos(ptri->angle[2])*(Eva-Evc)/ptri->edge_len[1] - Evb)
759 /(ptri->edge_len[0]*sin(ptri->angle[2]))/e;
760 }
761 mt->mapping(&pzone->danode[C],&aux[C],0);
762 AutoDScalar munBC = mt->mob->ElecMob(pc,nc,TD,adtl::fmax(0,Epnb),fabs(Etnb),TD);
763 AutoDScalar mupBC = mt->mob->HoleMob(pc,nc,TD,adtl::fmax(0,Eppb),fabs(Etpb),TD);
764
765 mt->mapping(&pzone->danode[A],&aux[A],0);
766 AutoDScalar munBA = mt->mob->ElecMob(pa,na,TD,adtl::fmax(0,Epnb),fabs(Etnb),TD);
767 AutoDScalar mupBA = mt->mob->HoleMob(pa,na,TD,adtl::fmax(0,Eppb),fabs(Etpb),TD);
768
769 mun = 0.5*(munBC+munBA);
770 mup = 0.5*(mupBC+mupBA);
771 }
772 else
773 {
774 mun = 0.5*(aux[C].mun+aux[A].mun);
775 mup = 0.5*(aux[C].mup+aux[A].mup);
776 }
777 F[0] = 0.5*(aux[C].eps+aux[A].eps)*(Va-Vc)/ptri->edge_len[1]*ptri->d[1];
778 F[1] = mun*Jnb*Jn_scale*ptri->d[1];
779 F[2] = - mup*Jpb*Jp_scale*ptri->d[1];
780 F[3] = - QNFactor*ptri->d[1]*bn*qV(qcc,qca,ptri->edge_len[1]);
781 F[4] = QPFactor*ptri->d[1]*bp*qV(qvc,qva,ptri->edge_len[1]);
782 for(int i=0;i<5;i++)
783 for(int j=0;j<5;j++)
784 {
785 J1.m[5*i+j] = F[i].getADValue(0+j);
786 J2.m[5*i+j] = F[i].getADValue(5+j);
787 J3.m[5*i+j] = F[i].getADValue(10+j);
788 }
789 MatSetValues(*jtmp,5,index3,5,index1,J1.m,ADD_VALUES);
790 MatSetValues(*jtmp,5,index3,5,index2,J2.m,ADD_VALUES);
791 MatSetValues(*jtmp,5,index3,5,index3,J3.m,ADD_VALUES);
792
793 F[0] = - 0.5*(aux[C].eps+aux[A].eps)*(Va-Vc)/ptri->edge_len[1]*ptri->d[1];
794 F[1] = - mun*Jnb*Jn_scale*ptri->d[1];
795 F[2] = mup*Jpb*Jp_scale*ptri->d[1];
796 F[3] = - QNFactor*ptri->d[1]*bn*qV(qca,qcc,ptri->edge_len[1]);
797 F[4] = QPFactor*ptri->d[1]*bp*qV(qva,qvc,ptri->edge_len[1]);
798 for(int i=0;i<5;i++)
799 for(int j=0;j<5;j++)
800 {
801 J1.m[5*i+j] = F[i].getADValue(0+j);
802 J2.m[5*i+j] = F[i].getADValue(5+j);
803 J3.m[5*i+j] = F[i].getADValue(10+j);
804 }
805 MatSetValues(*jtmp,5,index1,5,index1,J1.m,ADD_VALUES);
806 MatSetValues(*jtmp,5,index1,5,index2,J2.m,ADD_VALUES);
807 MatSetValues(*jtmp,5,index1,5,index3,J3.m,ADD_VALUES);
808
809 //---------------------------------------------------------------------------
810 Set_Mat5_zero(J1);
811 S = 0.25*ptri->d[2]*ptri->edge_len[2] + 0.25*ptri->d[1]*ptri->edge_len[1];
812 J1.m[6] = - Ra.getADValue(1)*S;
813 J1.m[7] = - Ra.getADValue(2)*S;
814 J1.m[11] = - Ra.getADValue(1)*S;
815 J1.m[12] = - Ra.getADValue(2)*S;
816 MatSetValues(*jtmp,5,index1,5,index1,J1.m,ADD_VALUES);
817
818 //---------------------------------------------------------------------------
819 Set_Mat5_zero(J2);
820 S = 0.25*ptri->d[0]*ptri->edge_len[0] + 0.25*ptri->d[2]*ptri->edge_len[2];
821 J2.m[6] = - Rb.getADValue(4)*S;
822 J2.m[7] = - Rb.getADValue(5)*S;
823 J2.m[11] = - Rb.getADValue(4)*S;
824 J2.m[12] = - Rb.getADValue(5)*S;
825 MatSetValues(*jtmp,5,index2,5,index2,J2.m,ADD_VALUES);
826
827 //---------------------------------------------------------------------------
828 Set_Mat5_zero(J3);
829 S = 0.25*ptri->d[1]*ptri->edge_len[1] + 0.25*ptri->d[0]*ptri->edge_len[0];
830 J3.m[6] = - Rc.getADValue(7)*S;
831 J3.m[7] = - Rc.getADValue(8)*S;
832 J3.m[11] = - Rc.getADValue(7)*S;
833 J3.m[12] = - Rc.getADValue(8)*S;
834 MatSetValues(*jtmp,5,index3,5,index3,J3.m,ADD_VALUES);
835 }
836
837 //-----------------------------------------------------------------------------
838 // boundaries
839 //-----------------------------------------------------------------------------
840
841
842 void SMCZone::F1Q_qddm_inner(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs)
843 {
844 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
845 PetscScalar e = mt->e;
846 PetscScalar hbar = mt->hbar;
847 PetscScalar Vi = x[zofs[zone_index]+5*i+0]; //potential of node i
848 PetscScalar ni = x[zofs[zone_index]+5*i+1]; //electron density of node i
849 PetscScalar pi = x[zofs[zone_index]+5*i+2]; //hole density of node i
850 PetscScalar Eqc = x[zofs[zone_index]+5*i+3];
851 PetscScalar Eqv = x[zofs[zone_index]+5*i+4];
852
853 f[zofs[zone_index]+5*i+0] = f[zofs[zone_index]+5*i+0]/pcell->area + mt->e*((pi-aux[i].Na)+(aux[i].Nd-ni));
854 f[zofs[zone_index]+5*i+1] = f[zofs[zone_index]+5*i+1]/pcell->area;
855 f[zofs[zone_index]+5*i+2] = f[zofs[zone_index]+5*i+2]/pcell->area;
856 f[zofs[zone_index]+5*i+3] = f[zofs[zone_index]+5*i+3]/pcell->area + (Eqc+Vi+aux[i].affinity);
857 f[zofs[zone_index]+5*i+4] = f[zofs[zone_index]+5*i+4]/pcell->area + (Eqv+Vi+aux[i].affinity+aux[i].Eg);
858
859 if(ODE_F.TimeDependent==true)
860 {
861 //second order
862 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
863 {
864 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
865 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
866 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
867 f[zofs[zone_index]+5*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
868 f[zofs[zone_index]+5*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
869 }
870 else //first order
871 {
872 f[zofs[zone_index]+5*i+1] += -(ni-fs[i].n)/ODE_F.dt;
873 f[zofs[zone_index]+5*i+2] += -(pi-fs[i].p)/ODE_F.dt;
874 }
875 }
876 }
877
878
879 void SMCZone::F1Q_qddm_ombc(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
880 {
881 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
882 int z = zone_index;
883 int equ_num = 5;
884 int size = pzone->davcell.size();
885 PetscScalar e = mt->e;
886 PetscScalar hbar = mt->hbar;
887 PetscScalar Na = aux[i].Na;
888 PetscScalar Nd = aux[i].Nd;
889 PetscScalar Vi = x[zofs[zone_index]+5*i+0]; //potential of node i
890 PetscScalar Eqc = x[zofs[zone_index]+5*i+3];
891 PetscScalar Eqv = x[zofs[zone_index]+5*i+4];
892 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
893 PetscScalar nie = mt->band->nie(fs[i].T);
894 PetscScalar Nc = mt->band->Nc(fs[i].T);
895 PetscScalar electron_density,hole_density;
896
897 int om_equ;
898 for(int j=0;j<electrode.size();j++)
899 {
900 if(electrode[j]==pcell->bc_index-1) {om_equ=j;break;}
901 }
902 f[zofs[z]+5*i+0] = Vi- mt->kb*fs[i].T/e*asinh((Nd-Na)/(2*nie)) + aux[i].affinity + mt->kb*fs[i].T/e*log(Nc/nie)
903 -x[zofs[z]+equ_num*size+om_equ];
904
905 if(Na>Nd) //p-type
906 {
907 hole_density = (-(Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
908 electron_density = nie*nie/hole_density;
909 }
910 else //n-type
911 {
912 electron_density = ((Nd-Na)+sqrt((Nd-Na)*(Nd-Na)+4*nie*nie))/2.0;
913 hole_density = nie*nie/electron_density;
914 }
915 f[zofs[z]+5*i+1] = x[zofs[z]+5*i+1] - electron_density; //electron density
916 f[zofs[z]+5*i+2] = x[zofs[z]+5*i+2] - hole_density; //hole density
917 f[zofs[z]+5*i+3] = f[zofs[z]+5*i+3]/pcell->area + (Eqc+Vi+aux[i].affinity);
918 f[zofs[z]+5*i+4] = f[zofs[z]+5*i+4]/pcell->area + (Eqv+Vi+aux[i].affinity+aux[i].Eg);
919 //f[zofs[z]+5*i+3] = x[zofs[z]+5*i+3] + x[zofs[z]+5*i+0]*e + aux[i].affinity*e;
920 //f[zofs[z]+5*i+4] = x[zofs[z]+5*i+4] + x[zofs[z]+5*i+0]*e + aux[i].affinity*e + aux[i].Eg;
921
922 }
923
924
925 void SMCZone::F1Q_qddm_stkbc(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int>& zofs,DABC &bc)
926 {
927 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
928 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
929 int z = zone_index;
930 int equ_num = 5;
931 int size = pzone->davcell.size();
932 PetscScalar Vi = x[zofs[zone_index]+5*i+0]; //potential of node i
933 PetscScalar ni = x[zofs[zone_index]+5*i+1]; //electron density of node i
934 PetscScalar pi = x[zofs[zone_index]+5*i+2]; //hole density of node i
935 PetscScalar Eqc = x[zofs[zone_index]+5*i+3];
936 PetscScalar Eqv = x[zofs[zone_index]+5*i+4];
937 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
938
939 int stk_equ;
940 for(int j=0;j<electrode.size();j++)
941 if(electrode[j]==pcell->bc_index-1) {stk_equ=j;break;}
942 //Schotty Barrier Lowerring
943 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
944 //Schottky current
945 PetscScalar Fn = mt->band->SchottyJsn(ni,fs[i].T,pbc->WorkFunction-aux[i].affinity-deltaVB)*0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
946 PetscScalar Fp = mt->band->SchottyJsp(pi,fs[i].T,pbc->WorkFunction-aux[i].affinity+deltaVB)*0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
947
948 f[zofs[z]+5*i+0] = x[zofs[z]+5*i+0] + pbc->WorkFunction - deltaVB - x[zofs[z]+equ_num*size+stk_equ];
949 f[zofs[z]+5*i+1] = (f[zofs[z]+5*i+1]+Fn)/pcell->area;
950 f[zofs[z]+5*i+2] = (f[zofs[z]+5*i+2]-Fp)/pcell->area;
951 f[zofs[z]+5*i+3] = f[zofs[z]+5*i+3]/pcell->area + (Eqc+Vi+aux[i].affinity);
952 f[zofs[z]+5*i+4] = f[zofs[z]+5*i+4]/pcell->area + (Eqv+Vi+aux[i].affinity+aux[i].Eg);
953
954 if(ODE_F.TimeDependent==true)
955 {
956 //second order
957 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
958 {
959 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
960 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
961 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
962 f[zofs[zone_index]+5*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
963 f[zofs[zone_index]+5*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
964 }
965 else //first order
966 {
967 f[zofs[zone_index]+5*i+1] += -(ni-fs[i].n)/ODE_F.dt;
968 f[zofs[zone_index]+5*i+2] += -(pi-fs[i].p)/ODE_F.dt;
969 }
970 }
971 }
972
973
974 void SMCZone::F1Q_qddm_insulator_gate(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
975 {
976 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
977 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
978 int equ_num = 5;
979 int size = pzone->davcell.size();
980 PetscScalar e = mt->e;
981 PetscScalar eV = mt->pscale->s_eV;
982 PetscScalar hbar = mt->hbar;
983 PetscScalar m0 = mt->me;
984 PetscScalar Vi = x[zofs[zone_index]+5*i+0]; //potential of node i
985 PetscScalar ni = x[zofs[zone_index]+5*i+1]; //electron density of node i
986 PetscScalar pi = x[zofs[zone_index]+5*i+2]; //hole density of node i
987 PetscScalar Eqc = x[zofs[zone_index]+5*i+3];
988 PetscScalar Eqv = x[zofs[zone_index]+5*i+4];
989 PetscScalar grad_P = 0, grad_qpn=0, grad_qpp=0;
990 int ins_equ;
991 for(int j=0;j<electrode.size();j++)
992 if(electrode[j]==pcell->bc_index-1)
993 {ins_equ=j;break;}
994 for(int j=0;j<pcell->nb_num;j++)
995 {
996 int nb = pcell->nb_array[j];
997 PetscScalar Vj = x[zofs[zone_index]+5*nb+0]; //potential of nb node
998
999 if(j==0||j==pcell->nb_num-1)
1000 {
1001 //the poisson's equation on third boundary type
1002 PetscScalar vgate = x[zofs[zone_index]+equ_num*size+ins_equ] - pbc->WorkFunction;
1003 PetscScalar q = mt->e*pbc->QF; //sigma is the surface change density
1004 PetscScalar Thick = pbc->Thick;
1005 PetscScalar eps_ox = mt->eps0*pbc->eps;
1006 PetscScalar r=q + eps_ox/Thick*vgate;
1007 PetscScalar s=eps_ox/Thick;
1008 grad_P += 0.5*pcell->ilen[j]*(r-0.25*s*(3*Vi+Vj));
1009 //the WBK approx of Si/SiO2 interface
1010 PetscScalar b_nox = hbar*hbar/(6*e*0.14*m0);
1011 PetscScalar b_pox = hbar*hbar/(6*e*1.0*m0);
1012 PetscScalar x_np = hbar/sqrt(2*0.4*m0*3.15*eV);
1013 PetscScalar x_pp = hbar/sqrt(2*0.4*m0*4.50*eV);
1014 grad_qpn += -b_nox/x_np*0.5*pcell->ilen[j];
1015 grad_qpp += -b_pox/x_pp*0.5*pcell->ilen[j];
1016 }
1017 }
1018
1019 f[zofs[zone_index]+5*i+0] = (f[zofs[zone_index]+5*i+0]+grad_P)/pcell->area
1020 + mt->e*((pi-aux[i].Na)+(aux[i].Nd-ni));
1021 f[zofs[zone_index]+5*i+1] = f[zofs[zone_index]+5*i+1]/pcell->area;
1022 f[zofs[zone_index]+5*i+2] = f[zofs[zone_index]+5*i+2]/pcell->area;
1023 f[zofs[zone_index]+5*i+3] = (f[zofs[zone_index]+5*i+3]+QNFactor*grad_qpn)/pcell->area + (Eqc+Vi+aux[i].affinity);
1024 f[zofs[zone_index]+5*i+4] = (f[zofs[zone_index]+5*i+4]-QPFactor*grad_qpp)/pcell->area + (Eqv+Vi+aux[i].affinity+aux[i].Eg);
1025
1026 if(ODE_F.TimeDependent==true)
1027 {
1028 //second order
1029 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1030 {
1031 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1032 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1033 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1034 f[zofs[zone_index]+5*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
1035 f[zofs[zone_index]+5*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
1036 }
1037 else //first order
1038 {
1039 f[zofs[zone_index]+5*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1040 f[zofs[zone_index]+5*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1041 }
1042 }
1043 }
1044
1045
1046 void SMCZone::F1Q_qddm_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1047 ISZone *pz, int n)
1048 {
1049 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1050 InsulatorInterfaceBC *pbc = dynamic_cast<InsulatorInterfaceBC * >(bc.Get_pointer(pcell->bc_index-1));
1051 PetscScalar e = mt->e;
1052 PetscScalar eV = mt->pscale->s_eV;
1053 PetscScalar hbar = mt->hbar;
1054 PetscScalar m0 = mt->me;
1055
1056 PetscScalar Vi = x[zofs[zone_index]+5*i+0]; //potential of node i
1057 PetscScalar ni = x[zofs[zone_index]+5*i+1]; //electron density of node i
1058 PetscScalar pi = x[zofs[zone_index]+5*i+2]; //hole density of node i
1059 PetscScalar Eqc = x[zofs[zone_index]+5*i+3];
1060 PetscScalar Eqv = x[zofs[zone_index]+5*i+4];
1061 PetscScalar Na = aux[i].Na;
1062 PetscScalar Nd = aux[i].Nd;
1063 PetscScalar L = 0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
1064
1065 PetscScalar grad_P = 0, grad_qpn=0, grad_qpp=0;
1066 for(int j=0;j<pcell->nb_num;j++)
1067 {
1068 int nb = pcell->nb_array[j];
1069 PetscScalar Vj = x[zofs[zone_index]+5*nb+0]; //potential of nb node
1070 grad_P += aux[i].eps*pcell->elen[j]/pcell->ilen[j]*(Vj-Vi);
1071 if(j==0||j==pcell->nb_num-1)
1072 {
1073 //the WBK approx of Si/SiO2 interface
1074 PetscScalar b_nox = hbar*hbar/(6*e*0.14*m0);
1075 PetscScalar b_pox = hbar*hbar/(6*e*1.0*m0);
1076 PetscScalar x_np = hbar/sqrt(2*0.4*m0*3.15*eV);
1077 PetscScalar x_pp = hbar/sqrt(2*0.4*m0*4.50*eV);
1078 grad_qpn += -b_nox/x_np*0.5*pcell->ilen[j];
1079 grad_qpp += -b_pox/x_pp*0.5*pcell->ilen[j];
1080 }
1081 }
1082 //poisson's equation at interface
1083 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1084 for(int j=0;j<ncell->nb_num;j++)
1085 {
1086 int nb = ncell->nb_array[j];
1087 PetscScalar Vj_n = x[zofs[pz->pzone->zone_index]+nb]; //potential of nb node
1088 grad_P += pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j]*(Vj_n-Vi);
1089 }
1090
1091 f[zofs[zone_index]+5*i+0] = grad_P + mt->e*((pi-ni)+(Nd-Na))*pcell->area + pbc->QF*L;
1092 f[zofs[zone_index]+5*i+1] = f[zofs[zone_index]+5*i+1]/pcell->area;
1093 f[zofs[zone_index]+5*i+2] = f[zofs[zone_index]+5*i+2]/pcell->area;
1094 f[zofs[zone_index]+5*i+3] = (f[zofs[zone_index]+5*i+3]+QNFactor*grad_qpn)/pcell->area + (Eqc+Vi+aux[i].affinity);
1095 f[zofs[zone_index]+5*i+4] = (f[zofs[zone_index]+5*i+4]-QPFactor*grad_qpp)/pcell->area + (Eqv+Vi+aux[i].affinity+aux[i].Eg);
1096
1097 if(ODE_F.TimeDependent==true)
1098 {
1099 //second order
1100 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1101 {
1102 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1103 PetscScalar Tn = (2-r)/(1-r)*ni-1.0/(r*(1-r))*fs[i].n+(1-r)/r*fs[i].n_last;
1104 PetscScalar Tp = (2-r)/(1-r)*pi-1.0/(r*(1-r))*fs[i].p+(1-r)/r*fs[i].p_last;
1105 f[zofs[zone_index]+5*i+1] += -Tn/(ODE_F.dt_last+ODE_F.dt);
1106 f[zofs[zone_index]+5*i+2] += -Tp/(ODE_F.dt_last+ODE_F.dt);
1107 }
1108 else //first order
1109 {
1110 f[zofs[zone_index]+5*i+1] += -(ni-fs[i].n)/ODE_F.dt;
1111 f[zofs[zone_index]+5*i+2] += -(pi-fs[i].p)/ODE_F.dt;
1112 }
1113 }
1114 }
1115
1116
1117 void SMCZone::F1Q_om_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1118 {
1119 int equ_num = 5;
1120 int size = pzone->davcell.size();
1121 int bc_index = electrode[i];
1122 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
1123 PetscScalar e = mt->e;
1124
1125 //calculate the total current of ohmic electrode
1126 PetscScalar current=0;
1127 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1128 {
1129 int node=bc[bc_index].psegment->node_array[j];
1130 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1131 PetscScalar Vi = x[zofs[zone_index]+5*node+0]; //potential of node i
1132 //conduction current
1133 current += DeviceDepth*(f[zofs[zone_index]+5*node+1]-f[zofs[zone_index]+5*node+2]);
1134 for(int k=0;k<pcell->nb_num;k++)
1135 {
1136 int nb = pcell->nb_array[k];
1137 PetscScalar Vj = x[zofs[zone_index]+5*nb+0]; //potential of nb node
1138 //displacement current
1139 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1140 }
1141 }
1142
1143 //if the electrode connect to another electrode
1144 if(pbc->inner_connect!=-1)
1145 {
1146 int connect_to = pbc->inner_connect;
1147 int connect_zone = pbc->connect_zone;
1148 int connect_elec=-1;
1149 SMCZone * pz = dynamic_cast <SMCZone * > (pbc->pzonedata);
1150 for(int j=0;j<pz->electrode.size();j++)
1151 {
1152 if(bc[pz->electrode[j]].BCType==OhmicContact)
1153 {
1154 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pz->electrode[j]));
1155 if(om_bc->inner_connect==bc_index) connect_elec=j;
1156 }
1157 }
1158 int connect_index=zofs[connect_zone]+equ_num*pz->pzone->davcell.size()+connect_elec;
1159
1160 if(bc_index < connect_to)
1161 {
1162 f[zofs[zone_index]+equ_num*size+i]+=x[zofs[zone_index]+equ_num*size+i]; //voltage equation
1163 f[connect_index]+=current; //current equation
1164 }
1165 else
1166 {
1167 f[zofs[zone_index]+equ_num*size+i]+=current; //current equation
1168 f[connect_index]+=-x[zofs[zone_index]+equ_num*size+i];//voltage equation
1169 }
1170 pbc->Set_Current_new(current);
1171 }
1172 //the electrode has an external voltage circult
1173 else if(pbc->electrode_type==VoltageBC)
1174 {
1175 PetscScalar V = pbc->Vapp;
1176 PetscScalar R = pbc->R;
1177 PetscScalar C = pbc->C;
1178 PetscScalar L = pbc->L;
1179 PetscScalar In = pbc->current;
1180 PetscScalar Icn = pbc->cap_current;
1181 PetscScalar Pn = pbc->potential;
1182 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1183 #ifdef _EXP_R_
1184 if(R<1.0)
1185 f[zofs[zone_index]+equ_num*size+i]=(L/ODE_F.dt+R)*current + (P-V) + (L/ODE_F.dt+R)*C/ODE_F.dt*P
1186 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn - L/ODE_F.dt*(In+Icn);
1187 else //when R is large, more like a current cource (R->+inf)
1188 f[zofs[zone_index]+equ_num*size+i]=current +(P-V)/(L/ODE_F.dt+R) + C/ODE_F.dt*P
1189 -C/ODE_F.dt*Pn - L/ODE_F.dt*(In+Icn)/(L/ODE_F.dt+R);
1190 #else
1191 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
1192 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1193 #endif
1194 pbc->Set_Current_new(current);
1195 }
1196 //the electrode has an external current circult
1197 else if(pbc->electrode_type==CurrentBC)
1198 {
1199 PetscScalar I = pbc->Iapp;
1200 PetscScalar C = pbc->C;
1201 f[zofs[zone_index]+equ_num*size+i]=current - I;
1202 pbc->Set_Current_new(I);
1203 }
1204 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1205 }
1206
1207
1208 void SMCZone::F1Q_stk_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1209 {
1210 int equ_num = 5;
1211 int size = pzone->davcell.size();
1212 int bc_index = electrode[i];
1213 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
1214
1215 PetscScalar current=0;
1216 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1217 {
1218 int node = bc[bc_index].psegment->node_array[j];
1219 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1220 PetscScalar Vi = x[zofs[zone_index]+5*node+0]; //potential of node i
1221 PetscScalar ni = x[zofs[zone_index]+5*node+1]; //electron density of node i
1222 PetscScalar pi = x[zofs[zone_index]+5*node+2]; //hole density of node i
1223 mt->mapping(&pzone->danode[node],&aux[node],ODE_F.clock);
1224 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
1225 PetscScalar Jsn = mt->band->SchottyJsn(ni,fs[node].T,pbc->WorkFunction-aux[node].affinity-deltaVB);
1226 PetscScalar Jsp = mt->band->SchottyJsp(pi,fs[node].T,pbc->WorkFunction-aux[node].affinity+deltaVB);
1227 current += -Jsn*0.5*pcell->ilen[0]*DeviceDepth;
1228 current += -Jsp*0.5*pcell->ilen[0]*DeviceDepth;
1229 current += -Jsn*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
1230 current += -Jsp*0.5*pcell->ilen[pcell->nb_num-1]*DeviceDepth;
1231 for(int k=0;k<pcell->nb_num;k++)
1232 {
1233 int nb = pcell->nb_array[k];
1234 PetscScalar Vj = x[zofs[zone_index]+5*nb+0]; //potential of nb node
1235 //displacement current
1236 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1237 }
1238 }
1239 if(pbc->electrode_type==VoltageBC)
1240 {
1241 PetscScalar R = pbc->R;
1242 PetscScalar C = pbc->C;
1243 PetscScalar L = pbc->L;
1244 PetscScalar V = pbc->Vapp;
1245 PetscScalar In = pbc->current;
1246 PetscScalar Icn = pbc->cap_current;
1247 PetscScalar Pn = pbc->potential;
1248 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1249 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
1250 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1251 pbc->Set_Current_new(current);
1252 }
1253 else if(pbc->electrode_type==CurrentBC)
1254 {
1255 PetscScalar I = pbc->Iapp;
1256 PetscScalar C = pbc->C;
1257 f[zofs[zone_index]+equ_num*size+i]=current - I;
1258 pbc->Set_Current_new(I);
1259 }
1260 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1261 }
1262
1263
1264 void SMCZone::F1Q_ins_electrode(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1265 {
1266 int equ_num = 5;
1267 int size = pzone->davcell.size();
1268 int bc_index = electrode[i];
1269 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
1270
1271 PetscScalar current=0;
1272 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1273 {
1274 int node = bc[bc_index].psegment->node_array[j];
1275 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1276 PetscScalar Vi = x[zofs[zone_index]+5*node+0]; //potential of node i
1277 for(int k=0;k<pcell->nb_num;k++)
1278 {
1279 int nb = pcell->nb_array[k];
1280 PetscScalar Vj = x[zofs[zone_index]+5*nb+0]; //potential of nb node
1281 //displacement current
1282 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
1283 }
1284 }
1285 if(pbc->electrode_type==VoltageBC)
1286 {
1287 PetscScalar R = pbc->R;
1288 PetscScalar C = pbc->C;
1289 PetscScalar L = pbc->L;
1290 PetscScalar V = pbc->Vapp;
1291 PetscScalar In = pbc->current;
1292 PetscScalar Icn = pbc->cap_current;
1293 PetscScalar Pn = pbc->potential;
1294 PetscScalar P = x[zofs[zone_index]+equ_num*size+i];
1295 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
1296 -(L/ODE_F.dt+R)*C/ODE_F.dt*Pn-L/ODE_F.dt*(In+Icn);
1297 pbc->Set_Current_new(current);
1298 }
1299 pbc->Set_Potential_new(x[zofs[zone_index]+equ_num*size+i]);
1300 }
1301
1302
1303
1304 void SMCZone::J1Q_qddm_inner(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs)
1305 {
1306 Mat5 A;
1307 PetscInt index[5],col[5];
1308 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1309 int z = zone_index;
1310
1311 PetscScalar Vi = x[zofs[z]+5*i+0]; //potential of node i
1312 PetscScalar ni = x[zofs[z]+5*i+1]; //electron density of node i
1313 PetscScalar pi = x[zofs[z]+5*i+2]; //hole density of node i
1314 PetscScalar area = pcell->area;
1315 PetscScalar e = mt->e;
1316
1317 //--------------------------------
1318 index[0] = zofs[z]+5*i+0;
1319 index[1] = zofs[z]+5*i+1;
1320 index[2] = zofs[z]+5*i+2;
1321 index[3] = zofs[z]+5*i+3;
1322 index[4] = zofs[z]+5*i+4;
1323
1324 //--------------------------------
1325 for(int j=0;j<pcell->nb_num;j++)
1326 {
1327 int nb = pcell->nb_array[j];
1328 col[0] = zofs[z]+5*nb+0;
1329 col[1] = zofs[z]+5*nb+1;
1330 col[2] = zofs[z]+5*nb+2;
1331 col[3] = zofs[z]+5*nb+3;
1332 col[4] = zofs[z]+5*nb+4;
1333 MatGetValues(*jtmp,5,index,5,col,A.m);
1334 A=A/area;
1335 MatSetValues(*jac,5,index,5,col,A.m,INSERT_VALUES);
1336 }
1337
1338 MatGetValues(*jtmp,5,index,5,index,A.m);
1339 A=A/area;
1340
1341 A.m[1] += -mt->e; //dfun(0)/dn(i)
1342 A.m[2] += mt->e; //dfun(0)/dp(i)
1343
1344 A.m[15] += 1;
1345 A.m[18] += 1;
1346
1347 A.m[20] += 1;
1348 A.m[24] += 1;
1349
1350 if(ODE_F.TimeDependent==true)
1351 {
1352 //second order
1353 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1354 {
1355 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1356 A.m[6] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1357 A.m[12] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1358 }
1359 else //first order
1360 {
1361 A.m[6] += -1/ODE_F.dt;
1362 A.m[12] += -1/ODE_F.dt;
1363 }
1364 }
1365
1366 MatSetValues(*jac,5,index,5,index,A.m,INSERT_VALUES);
1367 }
1368
1369
1370 void SMCZone::J1Q_qddm_ombc(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
1371 {
1372 Mat5 A;
1373 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1374 PetscScalar e = mt->e;
1375 PetscScalar area = pcell->area;
1376 PetscInt index[5],col[5];
1377 index[0] = zofs[zone_index]+5*i+0;
1378 index[1] = zofs[zone_index]+5*i+1;
1379 index[2] = zofs[zone_index]+5*i+2;
1380 index[3] = zofs[zone_index]+5*i+3;
1381 index[4] = zofs[zone_index]+5*i+4;
1382 //--------------------------------
1383 for(int j=0;j<pcell->nb_num;j++)
1384 {
1385 int nb = pcell->nb_array[j];
1386 col[0] = zofs[zone_index]+5*nb+0;
1387 col[1] = zofs[zone_index]+5*nb+1;
1388 col[2] = zofs[zone_index]+5*nb+2;
1389 col[3] = zofs[zone_index]+5*nb+3;
1390 col[4] = zofs[zone_index]+5*nb+4;
1391 MatGetValues(*jtmp,5,index,5,col,A.m);
1392 A=A/area;
1393 A.m[0]=0.0;A.m[1]=0.0;A.m[2]=0.0;A.m[3]=0.0;A.m[4]=0.0;
1394 A.m[5]=0.0;A.m[6]=0.0;A.m[7]=0.0;A.m[8]=0.0;A.m[9]=0.0;
1395 A.m[10]=0.0;A.m[11]=0.0;A.m[12]=0.0;A.m[13]=0.0;A.m[14]=0.0;
1396 MatSetValues(*jac,5,index,5,col,A.m,INSERT_VALUES);
1397 }
1398
1399 MatGetValues(*jtmp,5,index,5,index,A.m);
1400 A=A/area;
1401
1402 A.m[0] = 1.0;
1403 A.m[1] = 1.0;
1404 A.m[2] = 1.0;
1405
1406 A.m[5] = 0.0;
1407 A.m[6] = 1.0;
1408 A.m[7] = 0.0;
1409 A.m[8] = 0.0;
1410 A.m[9] = 0.0;
1411
1412 A.m[10] = 0.0;
1413 A.m[11] = 0.0;
1414 A.m[12] = 1.0;
1415 A.m[13] = 0.0;
1416 A.m[14] = 0.0;
1417
1418 A.m[15] += 1;
1419 A.m[18] += 1;
1420
1421 A.m[20] += 1;
1422 A.m[24] += 1;
1423
1424 MatSetValues(*jac,5,index,5,index,A.m,INSERT_VALUES);
1425 int om_equ;
1426 int equ_num = 5;
1427 int size = pzone->davcell.size();
1428 for(int j=0;j<electrode.size();j++)
1429 if(electrode[j]==pcell->bc_index-1) {om_equ=j;break;}
1430 MatSetValue(*jac,zofs[zone_index]+5*i,zofs[zone_index]+equ_num*size+om_equ,-1,INSERT_VALUES);
1431 }
1432
1433
1434 void SMCZone::J1Q_qddm_stkbc(int i,PetscScalar *x,Mat *jac,Mat *jtmp,ODE_Formula &ODE_F, vector<int> &zofs,DABC &bc)
1435 {
1436 Mat5 A;
1437 PetscInt index[5],col[5];
1438 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1439 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
1440 int z = zone_index;
1441 int equ_num = 5;
1442 int size = pzone->davcell.size();
1443 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
1444 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[i].eps,sqrt(aux[i].Ex*aux[i].Ex+aux[i].Ey*aux[i].Ey));
1445 PetscScalar Vi = x[zofs[z]+5*i+0]; //potential of node i
1446 PetscScalar ni = x[zofs[z]+5*i+1]; //electron density of node i
1447 PetscScalar pi = x[zofs[z]+5*i+2]; //hole density of node i
1448
1449 PetscScalar area = pcell->area;
1450 PetscScalar d_Fn_dni = 0;
1451 PetscScalar d_Fp_dpi = 0;
1452 //--------------------------------
1453 index[0] = zofs[z]+5*i+0;
1454 index[1] = zofs[z]+5*i+1;
1455 index[2] = zofs[z]+5*i+2;
1456 index[3] = zofs[z]+5*i+3;
1457 index[4] = zofs[z]+5*i+4;
1458 //--------------------------------
1459 for(int j=0;j<pcell->nb_num;j++)
1460 {
1461 int nb = pcell->nb_array[j];
1462 col[0] = zofs[z]+5*nb+0;
1463 col[1] = zofs[z]+5*nb+1;
1464 col[2] = zofs[z]+5*nb+2;
1465 col[3] = zofs[z]+5*nb+3;
1466 col[4] = zofs[z]+5*nb+4;
1467 MatGetValues(*jtmp,5,index,5,col,A.m);
1468 A.m[0] = 0.0;
1469 A=A/area;
1470 if(j==0||j==pcell->nb_num-1)
1471 {
1472 d_Fn_dni += mt->band->pdSchottyJsn_pdn(ni,fs[i].T,VB-deltaVB)*0.5*pcell->ilen[j];
1473 d_Fp_dpi += mt->band->pdSchottyJsp_pdp(pi,fs[i].T,VB+deltaVB)*0.5*pcell->ilen[j];
1474 }
1475 MatSetValues(*jac,5,index,5,col,A.m,INSERT_VALUES);
1476 }
1477
1478 MatGetValues(*jtmp,5,index,5,index,A.m);
1479 A=A/area;
1480 A.m[0] = 1.0; //dfun(0)/dP(i)
1481 A.m[1] = 0.0; //dfun(0)/dn(i)
1482 A.m[2] = 0.0; //dfun(0)/dp(i)
1483
1484 A.m[6] += d_Fn_dni/area; //dfun(1)/dn(i)
1485 A.m[12] += -d_Fp_dpi/area; //dfun(2)/dp(i)
1486
1487 A.m[15] += 1;
1488 A.m[18] += 1;
1489
1490 A.m[20] += 1;
1491 A.m[24] += 1;
1492
1493 if(ODE_F.TimeDependent==true)
1494 {
1495 //second order
1496 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1497 {
1498 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1499 A.m[6] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1500 A.m[12] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1501 }
1502 else //first order
1503 {
1504 A.m[6] += -1/ODE_F.dt;
1505 A.m[12] += -1/ODE_F.dt;
1506 }
1507 }
1508 MatSetValues(*jac,5,index,5,index,A.m,INSERT_VALUES);
1509
1510 int stk_equ;
1511 for(int j=0;j<electrode.size();j++)
1512 if(electrode[j]==pcell->bc_index-1)
1513 {stk_equ=j;break;}
1514 MatSetValue(*jac,zofs[zone_index]+5*i,zofs[zone_index]+equ_num*size+stk_equ,-1,INSERT_VALUES);
1515 }
1516
1517
1518 void SMCZone::J1Q_qddm_insulator_gate(int i,PetscScalar *x, Mat *jac, Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc)
1519 {
1520 Mat5 A;
1521 PetscInt index[5],col[5];
1522 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1523 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
1524 int z = zone_index;
1525 int equ_num = 5;
1526 int size = pzone->davcell.size();
1527
1528 PetscScalar ni = x[zofs[z]+5*i+1]; //electron density of node i
1529 PetscScalar pi = x[zofs[z]+5*i+2]; //hole density of node i
1530 PetscScalar e = mt->e;
1531 PetscScalar area = pcell->area;
1532 PetscScalar d_grad_P_dVi = 0;
1533 PetscScalar d_grad_P_dVapp = 0;
1534 int ins_equ;
1535 for(int j=0;j<electrode.size();j++)
1536 if(electrode[j]==pcell->bc_index-1)
1537 {ins_equ=j;break;}
1538 //--------------------------------
1539 index[0] = zofs[z]+5*i+0;
1540 index[1] = zofs[z]+5*i+1;
1541 index[2] = zofs[z]+5*i+2;
1542 index[3] = zofs[z]+5*i+3;
1543 index[4] = zofs[z]+5*i+4;
1544 //--------------------------------
1545 for(int j=0;j<pcell->nb_num;j++)
1546 {
1547 int nb = pcell->nb_array[j];
1548 col[0] = zofs[z]+5*nb+0;
1549 col[1] = zofs[z]+5*nb+1;
1550 col[2] = zofs[z]+5*nb+2;
1551 col[3] = zofs[z]+5*nb+3;
1552 col[4] = zofs[z]+5*nb+4;
1553 MatGetValues(*jtmp,5,index,5,col,A.m);
1554 A=A/area;
1555
1556 if(j==0||j==pcell->nb_num-1)
1557 {
1558 PetscScalar Thick = pbc->Thick;
1559 PetscScalar eps_ox = mt->eps0*pbc->eps;
1560 PetscScalar s=eps_ox/Thick;
1561 d_grad_P_dVi += -0.5*pcell->ilen[j]*0.25*s*3;
1562 d_grad_P_dVapp += 0.5*pcell->ilen[j]*s/area;
1563 A.m[0]+= -0.5*pcell->ilen[j]*0.25*s/area;
1564 }
1565 MatSetValues(*jac,5,index,5,col,A.m,INSERT_VALUES);
1566 }
1567
1568 MatGetValues(*jtmp,5,index,5,index,A.m);
1569 A=A/area;
1570 A.m[0] += d_grad_P_dVi/area; //dfun(0)/dP(i)
1571 A.m[1] += -mt->e; //dfun(0)/dn(i)
1572 A.m[2] += mt->e; //dfun(0)/dp(i)
1573
1574 A.m[15] += 1;
1575 A.m[18] += 1;
1576
1577 A.m[20] += 1;
1578 A.m[24] += 1;
1579
1580 if(ODE_F.TimeDependent==true)
1581 {
1582 //second order
1583 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1584 {
1585 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1586 A.m[6] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1587 A.m[12] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1588 }
1589 else //first order
1590 {
1591 A.m[6] += -1/ODE_F.dt;
1592 A.m[12] += -1/ODE_F.dt;
1593 }
1594 }
1595 MatSetValues(*jac,5,index,5,index,A.m,INSERT_VALUES);
1596 MatSetValue(*jac,zofs[z]+5*i,zofs[z]+equ_num*size+ins_equ,d_grad_P_dVapp,INSERT_VALUES);
1597 }
1598
1599
1600 void SMCZone::J1Q_qddm_interface(int i,PetscScalar *x,Mat *jac, Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc,
1601 ISZone *pz, int n)
1602 {
1603 Mat5 A;
1604 PetscInt index[5],col[5];
1605 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
1606 int z = zone_index;
1607 PetscScalar Vi = x[zofs[z]+5*i+0]; //potential of node i
1608 PetscScalar area = pcell->area;
1609 PetscScalar L = 0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
1610 PetscScalar d_grad_P_dVi = 0;
1611 //--------------------------------
1612 index[0] = zofs[z]+5*i+0;
1613 index[1] = zofs[z]+5*i+1;
1614 index[2] = zofs[z]+5*i+2;
1615 index[3] = zofs[z]+5*i+3;
1616 index[4] = zofs[z]+5*i+4;
1617 //--------------------------------
1618 for(int j=0;j<pcell->nb_num;j++)
1619 {
1620 int nb = pcell->nb_array[j];
1621 PetscScalar Vj = x[zofs[z]+5*nb+0]; //potential of nb node
1622 col[0] = zofs[z]+5*nb+0;
1623 col[1] = zofs[z]+5*nb+1;
1624 col[2] = zofs[z]+5*nb+2;
1625 col[3] = zofs[z]+5*nb+3;
1626 col[4] = zofs[z]+5*nb+4;
1627 //-------------------------------------
1628 d_grad_P_dVi += -aux[i].eps*pcell->elen[j]/pcell->ilen[j];
1629 MatGetValues(*jtmp,5,index,5,col,A.m);
1630 A=A/area;
1631 //-------------------------------------
1632 //df(x)i/dx(r)
1633 A.m[0] = aux[i].eps*pcell->elen[j]/pcell->ilen[j];
1634 A.m[1] = 0;
1635 A.m[2] = 0;
1636 MatSetValues(*jac,5,index,5,col,A.m,INSERT_VALUES);
1637 }
1638
1639 MatGetValues(*jtmp,5,index,5,index,A.m);
1640 A=A/area;
1641
1642 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
1643 int column = zofs[pz->pzone->zone_index]+n;
1644 for(int j=0;j<ncell->nb_num;j++)
1645 {
1646 int nb = ncell->nb_array[j];
1647 PetscScalar Vr_n = x[zofs[pz->pzone->zone_index]+nb]; //potential of nb node
1648 d_grad_P_dVi += -pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
1649 PetscScalar value = pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
1650 MatSetValue(*jac,zofs[z]+5*i+0,zofs[pz->pzone->zone_index]+nb,value,INSERT_VALUES);
1651 }
1652
1653 A.m[0] = d_grad_P_dVi; //dfun(0)/dP(i)
1654 A.m[1] = -mt->e*area; //dfun(0)/dn(i)
1655 A.m[2] = mt->e*area; //dfun(0)/dp(i)
1656
1657 A.m[15] += 1;
1658 A.m[18] += 1;
1659
1660 A.m[20] += 1;
1661 A.m[24] += 1;
1662
1663 if(ODE_F.TimeDependent==true)
1664 {
1665 //second order
1666 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
1667 {
1668 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
1669 A.m[6] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1670 A.m[12] += -(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
1671 }
1672 else //first order
1673 {
1674 A.m[6] += -1/ODE_F.dt;
1675 A.m[12] += -1/ODE_F.dt;
1676 }
1677 }
1678
1679 MatSetValues(*jac,5,index,5,index,A.m,INSERT_VALUES);
1680
1681 }
1682
1683
1684 void SMCZone::J1Q_om_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1685 {
1686 int equ_num = 5;
1687 int size = pzone->davcell.size();
1688 int bc_index = electrode[i];
1689 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
1690 PetscScalar A1[5],A2[5],J[5];
1691 PetscInt index[5],col[5];
1692 int matrix_row = zofs[zone_index]+equ_num*size+i;
1693 int connect_index;
1694 if(pbc->inner_connect!=-1)
1695 {
1696 int connect_to = pbc->inner_connect;
1697 int connect_zone = pbc->connect_zone;
1698 int connect_elec=-1;
1699 SMCZone * pz = dynamic_cast <SMCZone * > (pbc->pzonedata);
1700 for(int j=0;j<pz->electrode.size();j++)
1701 {
1702 if(bc[pz->electrode[j]].BCType==OhmicContact)
1703 {
1704 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pz->electrode[j]));
1705 if(om_bc->inner_connect==bc_index) connect_elec=j;
1706 }
1707 }
1708 connect_index=zofs[connect_zone]+equ_num*pz->pzone->davcell.size()+connect_elec;
1709 }
1710 PetscScalar R = pbc->R;
1711 PetscScalar C = pbc->C;
1712 PetscScalar L = pbc->L;
1713 PetscScalar e = mt->e;
1714 PetscScalar dEc_dV = -e;
1715 PetscScalar dEv_dV = -e;
1716
1717 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
1718 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
1719 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1720 {
1721 int node=bc[bc_index].psegment->node_array[j];
1722 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1723 PetscScalar dJdisp_dVi=0,dJdisp_dVr=0;
1724 index[0] = zofs[zone_index]+5*node+0;
1725 index[1] = zofs[zone_index]+5*node+1;
1726 index[2] = zofs[zone_index]+5*node+2;
1727 index[3] = zofs[zone_index]+5*node+3;
1728 index[4] = zofs[zone_index]+5*node+4;
1729 for(int k=0;k<pcell->nb_num;k++)
1730 {
1731 int nb = pcell->nb_array[k];
1732 col[0] = zofs[zone_index]+5*nb+0;
1733 col[1] = zofs[zone_index]+5*nb+1;
1734 col[2] = zofs[zone_index]+5*nb+2;
1735 col[3] = zofs[zone_index]+5*nb+3;
1736 col[4] = zofs[zone_index]+5*nb+4;
1737
1738 MatGetValues(*jtmp,1,&index[1],5,col,A1);
1739 MatGetValues(*jtmp,1,&index[2],5,col,A2);
1740
1741 //for displacement current
1742 dJdisp_dVi += aux[node].eps/pcell->ilen[k]/ODE_F.dt*pcell->elen[k];
1743 dJdisp_dVr = -aux[node].eps/pcell->ilen[k]/ODE_F.dt*pcell->elen[k];
1744
1745 //if the electrode connect to another electrode
1746 if(pbc->inner_connect!=-1)
1747 {
1748 int connect_to = pbc->inner_connect;
1749 if(bc_index < connect_to)
1750 {
1751 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1752 J[1]=(A1[1]-A2[1])*DeviceDepth;
1753 J[2]=(A1[2]-A2[2])*DeviceDepth;
1754 J[3]=(A1[3]-A2[3])*DeviceDepth;
1755 J[4]=(A1[4]-A2[4])*DeviceDepth;
1756 MatSetValues(*jac,1,&connect_index,5,col,J,ADD_VALUES);
1757 }
1758 else
1759 {
1760 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1761 J[1]=(A1[1]-A2[1])*DeviceDepth;
1762 J[2]=(A1[2]-A2[2])*DeviceDepth;
1763 J[3]=(A1[3]-A2[3])*DeviceDepth;
1764 J[4]=(A1[4]-A2[4])*DeviceDepth;
1765 MatSetValues(*jac,1,&matrix_row,5,col,J,ADD_VALUES);
1766 }
1767 }
1768 //the electrode has an external voltage circuit
1769 else if(pbc->electrode_type==VoltageBC)
1770 {
1771 #ifdef _EXP_R_
1772 if(R<1.0)
1773 {
1774 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth*(L/ODE_F.dt+R);
1775 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
1776 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
1777 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
1778 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
1779 }
1780 else
1781 {
1782 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1783 J[1]=(A1[1]-A2[1])*DeviceDepth;
1784 J[2]=(A1[2]-A2[2])*DeviceDepth;
1785 J[3]=(A1[3]-A2[3])*DeviceDepth;
1786 J[4]=(A1[4]-A2[4])*DeviceDepth;
1787 }
1788 #else
1789 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth*(L/ODE_F.dt+R);
1790 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
1791 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
1792 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
1793 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
1794 #endif
1795 MatSetValues(*jac,1,&matrix_row,5,col,J,ADD_VALUES);
1796 }
1797 //the electrode has an external current circuit
1798 else if(pbc->electrode_type==CurrentBC)
1799 {
1800 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1801 J[1]=(A1[1]-A2[1])*DeviceDepth;
1802 J[2]=(A1[2]-A2[2])*DeviceDepth;
1803 J[3]=(A1[3]-A2[3])*DeviceDepth;
1804 J[4]=(A1[4]-A2[4])*DeviceDepth;
1805 MatSetValues(*jac,1,&matrix_row,5,col,J,ADD_VALUES);
1806 }
1807 }
1808
1809 MatGetValues(*jtmp,1,&index[1],5,index,A1);
1810 MatGetValues(*jtmp,1,&index[2],5,index,A2);
1811
1812 if(pbc->inner_connect!=-1)
1813 {
1814 int connect_to = pbc->inner_connect;
1815 if(bc_index < connect_to)
1816 {
1817 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1818 J[1]=(A1[1]-A2[1])*DeviceDepth;
1819 J[2]=(A1[2]-A2[2])*DeviceDepth;
1820 J[3]=(A1[3]-A2[3])*DeviceDepth;
1821 J[4]=(A1[4]-A2[4])*DeviceDepth;
1822 MatSetValues(*jac,1,&connect_index,5,index,J,ADD_VALUES);
1823 }
1824 else
1825 {
1826 J[0]=(A1[0]-A2[0]+dJdisp_dVr)*DeviceDepth;
1827 J[1]=(A1[1]-A2[1])*DeviceDepth;
1828 J[2]=(A1[2]-A2[2])*DeviceDepth;
1829 J[3]=(A1[3]-A2[3])*DeviceDepth;
1830 J[4]=(A1[4]-A2[4])*DeviceDepth;
1831 MatSetValues(*jac,1,&matrix_row,5,index,J,ADD_VALUES);
1832 }
1833 }
1834 else if(pbc->electrode_type==VoltageBC)
1835 {
1836 #ifdef _EXP_R_
1837 if(R<1.0)
1838 {
1839 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth*(L/ODE_F.dt+R);
1840 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
1841 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
1842 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
1843 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
1844 }
1845 else
1846 {
1847 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth;
1848 J[1]=(A1[1]-A2[1])*DeviceDepth;
1849 J[2]=(A1[2]-A2[2])*DeviceDepth;
1850 J[3]=(A1[3]-A2[3])*DeviceDepth;
1851 J[4]=(A1[4]-A2[4])*DeviceDepth;
1852 }
1853 #else
1854 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth*(L/ODE_F.dt+R);
1855 J[1]=(A1[1]-A2[1])*DeviceDepth*(L/ODE_F.dt+R);
1856 J[2]=(A1[2]-A2[2])*DeviceDepth*(L/ODE_F.dt+R);
1857 J[3]=(A1[3]-A2[3])*DeviceDepth*(L/ODE_F.dt+R);
1858 J[4]=(A1[4]-A2[4])*DeviceDepth*(L/ODE_F.dt+R);
1859 #endif
1860 MatSetValues(*jac,1,&matrix_row,5,index,J,ADD_VALUES);
1861 }
1862 else if(pbc->electrode_type==CurrentBC)
1863 {
1864 J[0]=(A1[0]-A2[0]+dJdisp_dVi)*DeviceDepth;
1865 J[1]=(A1[1]-A2[1])*DeviceDepth;
1866 J[2]=(A1[2]-A2[2])*DeviceDepth;
1867 J[3]=(A1[3]-A2[3])*DeviceDepth;
1868 J[4]=(A1[4]-A2[4])*DeviceDepth;
1869 MatSetValues(*jac,1,&matrix_row,5,index,J,ADD_VALUES);
1870 }
1871
1872 }
1873
1874 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
1875 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
1876
1877 if(pbc->inner_connect!=-1)
1878 {
1879 int connect_to = pbc->inner_connect;
1880 if(bc_index < connect_to)
1881 {
1882 MatSetValue(*jac,matrix_row,matrix_row,1.0,ADD_VALUES);
1883 }
1884 else
1885 {
1886 MatSetValue(*jac,connect_index,matrix_row,-1.0,ADD_VALUES);
1887 }
1888 }
1889 else if(pbc->electrode_type==VoltageBC)
1890 {
1891 #ifdef _EXP_R_
1892 if(R<1.0)
1893 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
1894 else
1895 MatSetValue(*jac,matrix_row,matrix_row,1/(L/ODE_F.dt+R)+C/ODE_F.dt,INSERT_VALUES); //dJ/dP
1896 #else
1897 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
1898 #endif
1899 }
1900 }
1901
1902
1903 void SMCZone::J1Q_stk_electrode(int i,PetscScalar *x,Mat *jac,Mat *jtmp, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1904 {
1905 int equ_num = 5;
1906 int size = pzone->davcell.size();
1907 int bc_index = electrode[i];
1908 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
1909 int matrix_row = zofs[zone_index]+equ_num*size+i;
1910 PetscScalar R = pbc->R;
1911 PetscScalar C = pbc->C;
1912 PetscScalar L = pbc->L;
1913 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1914 {
1915 int node = bc[bc_index].psegment->node_array[j];
1916 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1917 PetscScalar VB = pbc->WorkFunction-aux[node].affinity;
1918 PetscScalar deltaVB=mt->band->SchottyBarrierLowerring(aux[node].eps,sqrt(aux[node].Ex*aux[node].Ex+aux[node].Ey*aux[node].Ey));
1919 PetscScalar ni = x[zofs[zone_index]+5*node+1]; //electron density of node i
1920 PetscScalar pi = x[zofs[zone_index]+5*node+2]; //hole density of node i
1921 PetscScalar dJ_dni = -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[0]
1922 -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB-deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
1923 PetscScalar dJ_dpi = -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[0]
1924 -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB+deltaVB)*0.5*pcell->ilen[pcell->nb_num-1];
1925
1926 //for displacement current
1927 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
1928 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
1929 for(int k=0;k<pcell->nb_num;k++)
1930 {
1931 int nb = pcell->nb_array[k];
1932 PetscScalar Vj = x[zofs[zone_index]+5*nb+0]; //potential of nb node
1933 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1934 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1935 if(pbc->electrode_type==VoltageBC)
1936 {
1937 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+0,dI_dVi*(L/ODE_F.dt+R),ADD_VALUES);
1938 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*nb+0,dI_dVr*(L/ODE_F.dt+R),ADD_VALUES);
1939 }
1940 else if(pbc->electrode_type==CurrentBC)
1941 {
1942 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+0,dI_dVi,ADD_VALUES);
1943 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*nb+0,dI_dVr,ADD_VALUES);
1944 }
1945 }
1946 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
1947 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
1948
1949 if(pbc->electrode_type==VoltageBC)
1950 {
1951 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+1,DeviceDepth*dJ_dni*(L/ODE_F.dt+R),INSERT_VALUES);
1952 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+2,DeviceDepth*dJ_dpi*(L/ODE_F.dt+R),INSERT_VALUES);
1953 }
1954 else if(pbc->electrode_type==CurrentBC)
1955 {
1956 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+1,DeviceDepth*dJ_dni,INSERT_VALUES);
1957 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+2,DeviceDepth*dJ_dpi,INSERT_VALUES);
1958 }
1959 }
1960 if(pbc->electrode_type==VoltageBC)
1961 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
1962 }
1963
1964
1965 void SMCZone::J1Q_ins_electrode(int i,PetscScalar *x,Mat *jac, Mat *jtmp,ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
1966 {
1967 int equ_num = 5;
1968 int size = pzone->davcell.size();
1969 int bc_index = electrode[i];
1970 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
1971 int matrix_row = zofs[zone_index]+equ_num*size+i;
1972 PetscScalar R = pbc->R;
1973 PetscScalar C = pbc->C;
1974 PetscScalar L = pbc->L;
1975
1976 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
1977 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
1978 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1979 {
1980 int node = bc[bc_index].psegment->node_array[j];
1981 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1982 //for displacement current
1983 for(int k=0;k<pcell->nb_num;k++)
1984 {
1985 int nb = pcell->nb_array[k];
1986 PetscScalar dI_dVi = DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1987 PetscScalar dI_dVr = -DeviceDepth*pcell->elen[k]*aux[node].eps/pcell->ilen[k]/ODE_F.dt;
1988 if(pbc->electrode_type==VoltageBC)
1989 {
1990 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*node+0,dI_dVi*(L/ODE_F.dt+R),ADD_VALUES);
1991 MatSetValue(*jac,matrix_row,zofs[zone_index]+5*nb+0,dI_dVr*(L/ODE_F.dt+R),ADD_VALUES);
1992 }
1993 }
1994 }
1995 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
1996 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
1997 if(pbc->electrode_type==VoltageBC)
1998 MatSetValue(*jac,matrix_row,matrix_row,1+(L/ODE_F.dt+R)*C/ODE_F.dt,INSERT_VALUES); //dJ/dP
1999
2000 }
2001
2002 void SMCZone::F1Q_efield_update(PetscScalar *x,vector<int> & zofs, DABC &bc, vector<BZoneData *>zonedata)
2003 {
2004 PetscScalar P_x=0,P_y=0,w=0;
2005 //calculate Ex Ey with least-squares gradient construction
2006 for(int i=0; i<pzone->davcell.size();i++)
2007 {
2008 P_x=0,P_y=0,w=0;
2009 VoronoiCell *pcell = pzone->davcell.GetPointer(i);
2010 for(int k=0;k<pcell->nb_num;k++)
2011 {
2012 const VoronoiCell *ncell = pzone->davcell.GetPointer(pcell->nb_array[k]);
2013 PetscScalar dx = ncell->x - pcell->x;
2014 PetscScalar dy = ncell->y - pcell->y;
2015 PetscScalar dP = x[zofs[zone_index]+5*pcell->nb_array[k]+0] - x[zofs[zone_index]+5*i+0];
2016 w=1.0/sqrt(dx*dx+dy*dy);
2017 P_x+=w*w*dP*dx;
2018 P_y+=w*w*dP*dy;
2019 }
2020 if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
2021 {
2022 InsulatorInterfaceBC *pbc;
2023 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
2024 int n_zone = pbc->pinterface->Find_neighbor_zone_index(zone_index);
2025 int n_node = pbc->pinterface->Find_neighbor_node_index(zone_index,i);
2026 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
2027 const VoronoiCell* dcell = pz->pzone->davcell.GetPointer(n_node);
2028 for(int k=1;k<dcell->nb_num-1;k++)
2029 {
2030 const VoronoiCell *ncell = pz->pzone->davcell.GetPointer(dcell->nb_array[k]);
2031 PetscScalar dx = ncell->x - dcell->x;
2032 PetscScalar dy = ncell->y - dcell->y;
2033 PetscScalar dP = x[zofs[n_zone]+dcell->nb_array[k]]- x[zofs[n_zone]+n_node];
2034 w=1.0/sqrt(dx*dx+dy*dy);
2035 P_x+=w*w*dP*dx;
2036 P_y+=w*w*dP*dy;
2037 }
2038 }
2039 aux[i].Ex = -(pcell->sc*P_x-pcell->sb*P_y);
2040 aux[i].Ey = -(pcell->sa*P_y-pcell->sb*P_x);
2041 }
2042 }
2043
General-purpose Semiconductor Simulator