?  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 / ddm1ac / semiequ1ac.cc
blob96b1ec2f538edd44206c213caa6e32e07583f7c8
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: Nov 19, 2006 */
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 "vec3.h"
23 #include <complex>
24 using namespace std;
25
26
27 void SMCZone::AC1_ddm_inner(int i, PetscScalar omega, PetscScalar *x, Mat *jac, Vec *r, Mat *jtmp, vector<int> &zofs)
28 {
29 Mat3 A;
30 Vec3 b;
31 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
32 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
33 int z = zone_index;
34 int N;
35 MatGetSize(*jtmp,&N,PETSC_NULL);
36 PetscScalar Vi = x[zofs[z]+3*i+0]; //potential of node i
37 PetscScalar ni = x[zofs[z]+3*i+1]; //electron density of node i
38 PetscScalar pi = x[zofs[z]+3*i+2]; //hole density of node i
39 PetscScalar area = pcell->area;
40
41 //--------------------------------
42 IndexRe[0] = zofs[z]+3*i+0;
43 IndexRe[1] = zofs[z]+3*i+1;
44 IndexRe[2] = zofs[z]+3*i+2;
45 IndexIm[0] = N + zofs[z]+3*i+0;
46 IndexIm[1] = N + zofs[z]+3*i+1;
47 IndexIm[2] = N + zofs[z]+3*i+2;
48 //--------------------------------
49 for(int j=0;j<pcell->nb_num;j++)
50 {
51 int nb = pcell->nb_array[j];
52 ColRe[0] = zofs[z]+3*nb+0;
53 ColRe[1] = zofs[z]+3*nb+1;
54 ColRe[2] = zofs[z]+3*nb+2;
55 ColIm[0] = N + zofs[z]+3*nb+0;
56 ColIm[1] = N + zofs[z]+3*nb+1;
57 ColIm[2] = N + zofs[z]+3*nb+2;
58 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A.m);
59 A=A/area;
60 MatSetValues(*jac,3,IndexRe,3,ColRe,A.m,INSERT_VALUES);
61 MatSetValues(*jac,3,IndexIm,3,ColIm,A.m,INSERT_VALUES);
62 }
63
64 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A.m);
65 A=A/area;
66
67 A.m[1] += -mt->e; //dfun(0)/dn(i)
68 A.m[2] += mt->e; //dfun(0)/dp(i)
69
70 MatSetValues(*jac,3,IndexRe,3,IndexRe,A.m,INSERT_VALUES);
71 MatSetValues(*jac,3,IndexIm,3,IndexIm,A.m,INSERT_VALUES);
72 Set_Mat3_zero(A);
73 A.m[4] = -omega;
74 A.m[8] = -omega;
75 MatSetValues(*jac,3,IndexIm,3,IndexRe,A.m,INSERT_VALUES);
76 A = -1.0*A;
77 MatSetValues(*jac,3,IndexRe,3,IndexIm,A.m,INSERT_VALUES);
78 Set_Vec3_zero(b);
79 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
80 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
81 }
82
83
84 void SMCZone::AC1_ddm_ombc(int i, PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,vector<int> &zofs,DABC &bc)
85 {
86 Mat3 A;
87 Vec3 b;
88 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
89 PetscInt IndexRe[3],IndexIm[3];
90 int z = zone_index;
91 int N;
92 MatGetSize(*jtmp,&N,PETSC_NULL);
93 Set_Mat3_I(A);
94 IndexRe[0] = zofs[z]+3*i+0;
95 IndexRe[1] = zofs[z]+3*i+1;
96 IndexRe[2] = zofs[z]+3*i+2;
97 IndexIm[0] = N + zofs[z]+3*i+0;
98 IndexIm[1] = N + zofs[z]+3*i+1;
99 IndexIm[2] = N + zofs[z]+3*i+2;
100 MatSetValues(*jac,3,IndexRe,3,IndexRe,A.m,INSERT_VALUES);
101 MatSetValues(*jac,3,IndexIm,3,IndexIm,A.m,INSERT_VALUES);
102 int om_equ;
103 int equ_num = 3;
104 int size = pzone->davcell.size();
105 for(int j=0;j<electrode.size();j++)
106 if(electrode[j]==pcell->bc_index-1) {om_equ=j;break;}
107 MatSetValue(*jac,zofs[z]+3*i,zofs[z]+equ_num*size+om_equ,-1,INSERT_VALUES);
108 MatSetValue(*jac,N+zofs[z]+3*i,N+zofs[z]+equ_num*size+om_equ,-1,INSERT_VALUES);
109 Set_Vec3_zero(b);
110 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
111 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
112 }
113
114
115 void SMCZone::AC1_ddm_stkbc(int i, PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,vector<int> &zofs,DABC &bc)
116 {
117 Mat3 A;
118 Vec3 b;
119 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
120 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
121 SchottkyBC *pbc = dynamic_cast<SchottkyBC * >(bc.Get_pointer(pcell->bc_index-1));
122 PetscScalar VB = pbc->WorkFunction-aux[i].affinity;
123 int z = zone_index;
124 int N;
125 MatGetSize(*jtmp,&N,PETSC_NULL);
126 PetscScalar Vi = x[zofs[z]+3*i+0]; //potential of node i
127 PetscScalar ni = x[zofs[z]+3*i+1]; //electron density of node i
128 PetscScalar pi = x[zofs[z]+3*i+2]; //hole density of node i
129 PetscScalar area = pcell->area;
130 PetscScalar d_Fn_dni = 0;
131 PetscScalar d_Fp_dpi = 0;
132 int stk_equ;
133 int equ_num = 3;
134 int size = pzone->davcell.size();
135 for(int j=0;j<electrode.size();j++)
136 if(electrode[j]==pcell->bc_index-1)
137 {stk_equ=j;break;}
138 //--------------------------------
139 IndexRe[0] = zofs[z]+3*i+0;
140 IndexRe[1] = zofs[z]+3*i+1;
141 IndexRe[2] = zofs[z]+3*i+2;
142 IndexIm[0] = N + zofs[z]+3*i+0;
143 IndexIm[1] = N + zofs[z]+3*i+1;
144 IndexIm[2] = N + zofs[z]+3*i+2;
145 //--------------------------------
146 for(int j=0;j<pcell->nb_num;j++)
147 {
148 int nb = pcell->nb_array[j];
149 ColRe[0] = zofs[z]+3*nb+0;
150 ColRe[1] = zofs[z]+3*nb+1;
151 ColRe[2] = zofs[z]+3*nb+2;
152 ColIm[0] = N + zofs[z]+3*nb+0;
153 ColIm[1] = N + zofs[z]+3*nb+1;
154 ColIm[2] = N + zofs[z]+3*nb+2;
155 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A.m);
156 A=A/area;
157 A.m[0]=A.m[1]=A.m[2]=0.0;
158 MatSetValues(*jac,3,IndexRe,3,ColRe,A.m,INSERT_VALUES);
159 MatSetValues(*jac,3,IndexIm,3,ColIm,A.m,INSERT_VALUES);
160 if(j==0||j==pcell->nb_num-1)
161 {
162 d_Fn_dni += mt->band->pdSchottyJsn_pdn(ni,fs[i].T,VB)*0.5*pcell->ilen[j];
163 d_Fp_dpi += mt->band->pdSchottyJsp_pdp(pi,fs[i].T,VB)*0.5*pcell->ilen[j];
164 }
165 }
166
167 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A.m);
168 A=A/area;
169 A.m[0] = 1.0;
170 A.m[1] = 0.0; //dfun(0)/dn(i)
171 A.m[2] = 0.0; //dfun(0)/dp(i)
172
173
174 A.m[4] += d_Fn_dni/area; //dfun(1)/dn(i)
175 A.m[8] += - d_Fp_dpi/area; //dfun(2)/dp(i)
176
177 MatSetValues(*jac,3,IndexRe,3,IndexRe,A.m,INSERT_VALUES);
178 MatSetValues(*jac,3,IndexIm,3,IndexIm,A.m,INSERT_VALUES);
179 Set_Mat3_zero(A);
180 A.m[4] = -omega;
181 A.m[8] = -omega;
182 MatSetValues(*jac,3,IndexIm,3,IndexRe,A.m,INSERT_VALUES);
183 A = -1.0*A;
184 MatSetValues(*jac,3,IndexRe,3,IndexIm,A.m,INSERT_VALUES);
185
186 MatSetValue(*jac,zofs[z]+3*i,zofs[z]+equ_num*size+stk_equ,-1,INSERT_VALUES);
187 MatSetValue(*jac,N+zofs[z]+3*i,N+zofs[z]+equ_num*size+stk_equ,-1,INSERT_VALUES);
188 Set_Vec3_zero(b);
189 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
190 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
191 }
192
193
194 void SMCZone::AC1_ddm_insulator_gate(int i, PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,vector<int> &zofs,DABC &bc)
195 {
196 Mat3 A;
197 Vec3 b;
198 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
199 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
200 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
201 int z = zone_index;
202 int equ_num = 3;
203 int size = pzone->davcell.size();
204 int N;
205 MatGetSize(*jtmp,&N,PETSC_NULL);
206
207 PetscScalar ni = x[zofs[z]+3*i+1]; //electron density of node i
208 PetscScalar pi = x[zofs[z]+3*i+2]; //hole density of node i
209
210 PetscScalar area = pcell->area;
211 PetscScalar L = 0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
212 PetscScalar d_grad_P_dVi = 0;
213 PetscScalar d_grad_P_dVapp = 0;
214 int ins_equ;
215 for(int j=0;j<electrode.size();j++)
216 if(electrode[j]==pcell->bc_index-1)
217 {ins_equ=j;break;}
218 //--------------------------------
219 IndexRe[0] = zofs[z]+3*i+0;
220 IndexRe[1] = zofs[z]+3*i+1;
221 IndexRe[2] = zofs[z]+3*i+2;
222 IndexIm[0] = N + zofs[z]+3*i+0;
223 IndexIm[1] = N + zofs[z]+3*i+1;
224 IndexIm[2] = N + zofs[z]+3*i+2;
225 for(int j=0;j<pcell->nb_num;j++)
226 {
227 int nb = pcell->nb_array[j];
228 ColRe[0] = zofs[z]+3*nb+0;
229 ColRe[1] = zofs[z]+3*nb+1;
230 ColRe[2] = zofs[z]+3*nb+2;
231 ColIm[0] = N + zofs[z]+3*nb+0;
232 ColIm[1] = N + zofs[z]+3*nb+1;
233 ColIm[2] = N + zofs[z]+3*nb+2;
234 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A.m);
235 A=A/area;
236 if(j==0||j==pcell->nb_num-1)
237 {
238 PetscScalar Thick = pbc->Thick;
239 PetscScalar eps_ox = mt->eps0*pbc->eps;
240 PetscScalar s=eps_ox/Thick;
241 d_grad_P_dVi += -0.5*pcell->ilen[j]*0.25*s*3;
242 d_grad_P_dVapp += 0.5*pcell->ilen[j]*s/area;
243 A.m[0]+= -0.5*pcell->ilen[j]*0.25*s/area;
244 }
245 MatSetValues(*jac,3,IndexRe,3,ColRe,A.m,INSERT_VALUES);
246 MatSetValues(*jac,3,IndexIm,3,ColIm,A.m,INSERT_VALUES);
247 }
248 mt->mapping(&pzone->danode[i],&aux[i],0);
249
250 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A.m);
251 A=A/area;
252 A.m[0] += d_grad_P_dVi/area; //dfun(0)/dP(i)
253 A.m[1] += -mt->e; //dfun(0)/dn(i)
254 A.m[2] += mt->e; //dfun(0)/dp(i)
255
256 MatSetValues(*jac,3,IndexRe,3,IndexRe,A.m,INSERT_VALUES);
257 MatSetValues(*jac,3,IndexIm,3,IndexIm,A.m,INSERT_VALUES);
258
259 Set_Mat3_zero(A);
260 A.m[4] = -omega;
261 A.m[8] = -omega;
262 MatSetValues(*jac,3,IndexIm,3,IndexRe,A.m,INSERT_VALUES);
263 A = -1.0*A;
264 MatSetValues(*jac,3,IndexRe,3,IndexIm,A.m,INSERT_VALUES);
265
266 MatSetValue(*jac,zofs[z]+3*i,zofs[z]+equ_num*size+ins_equ,d_grad_P_dVapp,INSERT_VALUES);
267 MatSetValue(*jac,N+zofs[z]+3*i,N+zofs[z]+equ_num*size+ins_equ,d_grad_P_dVapp,INSERT_VALUES);
268
269 Set_Vec3_zero(b);
270 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
271 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
272 }
273
274
275 void SMCZone::AC1_ddm_interface(int i, PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,
276 vector<int> &zofs,DABC &bc,ISZone *pz, int n)
277 {
278 Mat3 A;
279 Vec3 b;
280 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
281 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
282 InsulatorContactBC *pbc = dynamic_cast<InsulatorContactBC * >(bc.Get_pointer(pcell->bc_index-1));
283 int z = zone_index;
284 int equ_num = 3;
285 int size = pzone->davcell.size();
286 int N;
287 MatGetSize(*jtmp,&N,PETSC_NULL);
288
289 PetscScalar ni = x[zofs[z]+3*i+1]; //electron density of node i
290 PetscScalar pi = x[zofs[z]+3*i+2]; //hole density of node i
291
292 PetscScalar area = pcell->area;
293 PetscScalar L = 0.5*(pcell->ilen[0]+pcell->ilen[pcell->nb_num-1]);
294 PetscScalar d_grad_P_dVi = 0;
295
296 //--------------------------------
297 IndexRe[0] = zofs[z]+3*i+0;
298 IndexRe[1] = zofs[z]+3*i+1;
299 IndexRe[2] = zofs[z]+3*i+2;
300 IndexIm[0] = N + zofs[z]+3*i+0;
301 IndexIm[1] = N + zofs[z]+3*i+1;
302 IndexIm[2] = N + zofs[z]+3*i+2;
303 for(int j=0;j<pcell->nb_num;j++)
304 {
305 int nb = pcell->nb_array[j];
306 ColRe[0] = zofs[z]+3*nb+0;
307 ColRe[1] = zofs[z]+3*nb+1;
308 ColRe[2] = zofs[z]+3*nb+2;
309 ColIm[0] = N + zofs[z]+3*nb+0;
310 ColIm[1] = N + zofs[z]+3*nb+1;
311 ColIm[2] = N + zofs[z]+3*nb+2;
312 d_grad_P_dVi += -aux[i].eps*pcell->elen[j]/pcell->ilen[j];
313 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A.m);
314 A=A/area;
315 A.m[0] = aux[i].eps*pcell->elen[j]/pcell->ilen[j];
316 A.m[1] = 0;
317 A.m[2] = 0;
318 MatSetValues(*jac,3,IndexRe,3,ColRe,A.m,INSERT_VALUES);
319 MatSetValues(*jac,3,IndexIm,3,ColIm,A.m,INSERT_VALUES);
320 }
321 mt->mapping(&pzone->danode[i],&aux[i],0);
322
323 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A.m);
324 A=A/area;
325 A.m[0] = d_grad_P_dVi; //dfun(0)/dP(i)
326 A.m[1] = -mt->e*area; //dfun(0)/dn(i)
327 A.m[2] = mt->e*area; //dfun(0)/dp(i)
328 MatSetValues(*jac,3,IndexRe,3,IndexRe,A.m,INSERT_VALUES);
329 MatSetValues(*jac,3,IndexIm,3,IndexIm,A.m,INSERT_VALUES);
330
331 Set_Mat3_zero(A);
332 A.m[4] = -omega;
333 A.m[8] = -omega;
334 MatSetValues(*jac,3,IndexIm,3,IndexRe,A.m,INSERT_VALUES);
335 A = -1.0*A;
336 MatSetValues(*jac,3,IndexRe,3,IndexIm,A.m,INSERT_VALUES);
337
338 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
339 int column = zofs[pz->pzone->zone_index]+n;
340 for(int j=0;j<ncell->nb_num;j++)
341 {
342 int nb = ncell->nb_array[j];
343 PetscScalar Vr_n = x[zofs[pz->pzone->zone_index]+nb]; //potential of nb node
344 d_grad_P_dVi += -pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
345 PetscScalar d_grad_P_dVj = pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
346 MatSetValue(*jac,zofs[z]+3*i+0,zofs[pz->pzone->zone_index]+nb,d_grad_P_dVj,INSERT_VALUES);
347 MatSetValue(*jac,N+zofs[z]+3*i+0,N+zofs[pz->pzone->zone_index]+nb,d_grad_P_dVj,INSERT_VALUES);
348 }
349 MatSetValue(*jac,zofs[z]+3*i+0,zofs[z]+3*i+0,d_grad_P_dVi,INSERT_VALUES);
350 MatSetValue(*jac,N+zofs[z]+3*i+0,N+zofs[z]+3*i+0,d_grad_P_dVi,INSERT_VALUES);
351
352 Set_Vec3_zero(b);
353 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
354 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
355 }
356
357
358 void SMCZone::AC1_ddm_homojunction(int i, PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,
359 vector<int> &zofs,DABC &bc,SMCZone *pz, int n)
360 {
361 Mat3 A1,A2,AJ;
362 Vec3 b;
363 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3],NbRe[3],NbIm[3];
364 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
365 int z = zone_index;
366 int N;
367 MatGetSize(*jtmp,&N,PETSC_NULL);
368 PetscScalar e = mt->e;
369 PetscScalar Vi = x[zofs[z]+3*i+0]; //potential of node i
370 PetscScalar ni = x[zofs[z]+3*i+1]; //electron density of node i
371 PetscScalar pi = x[zofs[z]+3*i+2]; //hole density of node i
372
373 PetscScalar d_grad_P_dVi = 0;
374 //--------------------------------
375 IndexRe[0] = zofs[z]+3*i+0;
376 IndexRe[1] = zofs[z]+3*i+1;
377 IndexRe[2] = zofs[z]+3*i+2;
378 IndexIm[0] = N + zofs[z]+3*i+0;
379 IndexIm[1] = N + zofs[z]+3*i+1;
380 IndexIm[2] = N + zofs[z]+3*i+2;
381 //--------------------------------
382
383 if(zone_index > pz->pzone->zone_index)
384 {
385 Set_Mat3_I(AJ);
386 MatSetValues(*jac,3,IndexRe,3,IndexRe,AJ.m,INSERT_VALUES);
387 MatSetValues(*jac,3,IndexIm,3,IndexIm,AJ.m,INSERT_VALUES);
388 AJ = -1.0*AJ;
389 NbRe[0] = zofs[pz->zone_index]+3*n+0;
390 NbRe[1] = zofs[pz->zone_index]+3*n+1;
391 NbRe[2] = zofs[pz->zone_index]+3*n+2;
392 NbIm[0] = N+zofs[pz->zone_index]+3*n+0;
393 NbIm[1] = N+zofs[pz->zone_index]+3*n+1;
394 NbIm[2] = N+zofs[pz->zone_index]+3*n+2;
395 MatSetValues(*jac,3,IndexRe,3,NbRe,AJ.m,INSERT_VALUES);
396 MatSetValues(*jac,3,IndexIm,3,NbIm,AJ.m,INSERT_VALUES);
397 return;
398 }
399
400 for(int j=0;j<pcell->nb_num;j++)
401 {
402 int nb = pcell->nb_array[j];
403 PetscScalar Vj = x[zofs[z]+3*nb+0]; //potential of nb node
404 //-------------------------------------
405 d_grad_P_dVi += -aux[i].eps*pcell->elen[j]/pcell->ilen[j];
406 ColRe[0] = zofs[z]+3*nb+0;
407 ColRe[1] = zofs[z]+3*nb+1;
408 ColRe[2] = zofs[z]+3*nb+2;
409 ColIm[0] = N+zofs[z]+3*nb+0;
410 ColIm[1] = N+zofs[z]+3*nb+1;
411 ColIm[2] = N+zofs[z]+3*nb+2;
412 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A1.m);
413 //-------------------------------------
414 //df(x)i/dx(r)
415 A1.m[0] = aux[i].eps*pcell->elen[j]/pcell->ilen[j]; //dfun(0)/dP(r)
416 A1.m[1] = 0; //dfun(0)/dn(r)
417 A1.m[2] = 0; //dfun(0)/dp(r)
418 MatSetValues(*jac,3,IndexRe,3,ColRe,A1.m,INSERT_VALUES);
419 MatSetValues(*jac,3,IndexIm,3,ColIm,A1.m,INSERT_VALUES);
420 }
421
422 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
423 NbRe[0] = zofs[pz->zone_index]+3*n+0;
424 NbRe[1] = zofs[pz->zone_index]+3*n+1;
425 NbRe[2] = zofs[pz->zone_index]+3*n+2;
426 NbIm[0] = N+zofs[pz->zone_index]+3*n+0;
427 NbIm[1] = N+zofs[pz->zone_index]+3*n+1;
428 NbIm[2] = N+zofs[pz->zone_index]+3*n+2;
429 for(int j=0;j<ncell->nb_num;j++)
430 {
431 int nb = ncell->nb_array[j];
432 PetscScalar Vj = x[zofs[pz->zone_index]+3*nb+0]; //potential of nb node
433 //-------------------------------------
434 d_grad_P_dVi += -pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
435 ColRe[0] = zofs[pz->zone_index]+3*nb+0;
436 ColRe[1] = zofs[pz->zone_index]+3*nb+1;
437 ColRe[2] = zofs[pz->zone_index]+3*nb+2;
438 ColIm[0] = N+zofs[pz->zone_index]+3*nb+0;
439 ColIm[1] = N+zofs[pz->zone_index]+3*nb+1;
440 ColIm[2] = N+zofs[pz->zone_index]+3*nb+2;
441 MatGetValues(*jtmp,3,NbRe,3,ColRe,A2.m);
442 //-------------------------------------
443 //df(x)i/dx(r)
444 A2.m[0] = pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j]; //dfun(0)/dP(r)
445 A2.m[1] = 0; //dfun(0)/dn(r)
446 A2.m[2] = 0; //dfun(0)/dp(r)
447 MatSetValues(*jac,3,IndexRe,3,ColRe,A2.m,INSERT_VALUES);
448 MatSetValues(*jac,3,IndexIm,3,ColIm,A2.m,INSERT_VALUES);
449 }
450
451 //fun(0) is the poisson's equation
452 //fun(1) is the continuous equation of electron
453 //fun(2) is the continuous equation of hole
454 PetscScalar area = pcell->area + ncell->area;
455 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A1.m);
456 MatGetValues(*jtmp,3,NbRe,3,NbRe,A2.m);
457 AJ=A1+A2;
458
459 AJ.m[0] = d_grad_P_dVi; //dfun(0)/dP(i)
460 AJ.m[1] = -e*area; //dfun(0)/dn(i)
461 AJ.m[2] = e*area; //dfun(0)/dp(i)
462
463 MatSetValues(*jac,3,IndexRe,3,IndexRe,AJ.m,INSERT_VALUES);
464 MatSetValues(*jac,3,IndexIm,3,IndexIm,AJ.m,INSERT_VALUES);
465 Set_Mat3_zero(AJ);
466 AJ.m[4] = -omega*area;
467 AJ.m[8] = -omega*area;
468 MatSetValues(*jac,3,IndexIm,3,IndexRe,AJ.m,INSERT_VALUES);
469 AJ = -1.0*AJ;
470 MatSetValues(*jac,3,IndexRe,3,IndexIm,AJ.m,INSERT_VALUES);
471 Set_Vec3_zero(b);
472 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
473 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
474 }
475
476
477 void SMCZone::AC1_ddm_heterojunction(int i, PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,
478 vector<int> &zofs,DABC &bc,SMCZone *pz, int n)
479 {
480 Mat3 A;
481 Vec3 b;
482 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
483 const VoronoiCell* pcell = pzone->davcell.GetPointer(i);
484 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
485 int z = zone_index;
486 int N;
487 MatGetSize(*jtmp,&N,PETSC_NULL);
488 PetscScalar k = mt->kb;
489 PetscScalar e = mt->e;
490 PetscScalar Vi = x[zofs[z]+3*i+0]; //potential of node i
491 PetscScalar ni = x[zofs[z]+3*i+1]; //electron density of node i
492 PetscScalar pi = x[zofs[z]+3*i+2]; //hole density of node i
493 PetscScalar Ti = fs[i].T;
494 PetscScalar Eci = -e*(Vi + aux[i].affinity + mt->kb*fs[i].T/e*log(aux[i].Nc));
495 PetscScalar Evi = -e*(Vi + aux[i].affinity + aux[i].Eg/e - mt->kb*fs[i].T/e*log(aux[i].Nv));
496 PetscScalar area = pcell->area;
497 PetscScalar d_grad_P_dVi = 0;
498 PetscScalar d_Fn_dni = 0;
499 PetscScalar d_Fp_dpi = 0;
500 //--------------------------------
501 IndexRe[0] = zofs[z]+3*i+0;
502 IndexRe[1] = zofs[z]+3*i+1;
503 IndexRe[2] = zofs[z]+3*i+2;
504 IndexIm[0] = N + zofs[z]+3*i+0;
505 IndexIm[1] = N + zofs[z]+3*i+1;
506 IndexIm[2] = N + zofs[z]+3*i+2;
507 //--------------------------------
508 for(int j=0;j<pcell->nb_num;j++)
509 {
510 int nb = pcell->nb_array[j];
511
512 PetscScalar Vj = x[zofs[z]+3*nb+0]; //potential of nb node
513 PetscScalar nj = x[zofs[z]+3*nb+1]; //electron density of nb node
514 PetscScalar pj = x[zofs[z]+3*nb+2]; //hole density of nb node
515 //-------------------------------------
516 d_grad_P_dVi += -aux[i].eps*pcell->elen[j]/pcell->ilen[j];
517
518 ColRe[0] = zofs[z]+3*nb+0;
519 ColRe[1] = zofs[z]+3*nb+1;
520 ColRe[2] = zofs[z]+3*nb+2;
521 ColIm[0] = N + zofs[z]+3*nb+0;
522 ColIm[1] = N + zofs[z]+3*nb+1;
523 ColIm[2] = N + zofs[z]+3*nb+2;
524 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A.m);
525 A=A/area;
526 //-------------------------------------
527 //df(x)i/dx(r)
528 if(pzone->zone_index < pz->pzone->zone_index)
529 {
530 A.m[0] = aux[i].eps*pcell->elen[j]/pcell->ilen[j];//dfun(0)/dP(r)
531 A.m[1] = 0; //dfun(0)/dn(r)
532 A.m[2] = 0; //dfun(0)/dp(r)
533 }
534 else
535 {
536 A.m[0] = 0; //dfun(0)/dP(r)
537 A.m[1] = 0; //dfun(0)/dn(r)
538 A.m[2] = 0; //dfun(0)/dp(r)
539 }
540
541 //Thermal emit current
542 if(j==0||j==pcell->nb_num-1)
543 {
544 PetscScalar Ecj = -e*(Vi + pz->aux[n].affinity + k*pz->fs[n].T/e*log(pz->aux[n].Nc));
545 PetscScalar Evj = -e*(Vi + pz->aux[n].affinity - k*pz->fs[n].T/e*log(pz->aux[n].Nv) + pz->aux[n].Eg/e);
546 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
547 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
548 if(Ecj > Eci)
549 {
550 PetscScalar pm = pz->mt->band->EffecElecMass(Ti)/mt->band->EffecElecMass(Ti);
551 d_Fn_dni += 2*e*pm*mt->band->ThermalVn(Ti)*exp(-(Ecj-Eci)/(k*Ti))*0.5*pcell->ilen[j]/pcell->area;
552 PetscScalar d_Fn_dnd = -2*e*pz->mt->band->ThermalVn(Ti)*0.5*pcell->ilen[j]/pcell->area;
553 MatSetValue(*jac,zofs[z]+3*i+1,zofs[pz->pzone->zone_index]+3*n+1,d_Fn_dnd,ADD_VALUES);
554 MatSetValue(*jac,N+zofs[z]+3*i+1,N+zofs[pz->pzone->zone_index]+3*n+1,d_Fn_dnd,ADD_VALUES);
555 }
556 else
557 {
558 PetscScalar pm = mt->band->EffecElecMass(Ti)/pz->mt->band->EffecElecMass(Ti);
559 d_Fn_dni += -2*e*mt->band->ThermalVn(Ti)*0.5*pcell->ilen[j]/pcell->area;
560 PetscScalar d_Fn_dnd = 2*e*pm*pz->mt->band->ThermalVn(Ti)*exp(-(Eci-Ecj)/(k*Ti))*0.5*pcell->ilen[j]/pcell->area;
561 MatSetValue(*jac,zofs[z]+3*i+1,zofs[pz->pzone->zone_index]+3*n+1,d_Fn_dnd,ADD_VALUES);
562 MatSetValue(*jac,N+zofs[z]+3*i+1,N+zofs[pz->pzone->zone_index]+3*n+1,d_Fn_dnd,ADD_VALUES);
563 }
564
565 if(Evj > Evi)
566 {
567 PetscScalar pm = pz->mt->band->EffecHoleMass(Ti)/mt->band->EffecHoleMass(Ti);
568 d_Fp_dpi += -2*e*pm*mt->band->ThermalVp(Ti)*exp(-(Evj-Evi)/(k*Ti))*0.5*pcell->ilen[j]/pcell->area;
569 PetscScalar d_Fp_dpd = 2*e*pz->mt->band->ThermalVp(Ti)*0.5*pcell->ilen[j]/pcell->area;
570 MatSetValue(*jac,zofs[z]+3*i+2,zofs[pz->pzone->zone_index]+3*n+2,-d_Fp_dpd,ADD_VALUES);
571 MatSetValue(*jac,N+zofs[z]+3*i+2,N+zofs[pz->pzone->zone_index]+3*n+2,-d_Fp_dpd,ADD_VALUES);
572 }
573 else
574 {
575 PetscScalar pm = mt->band->EffecHoleMass(Ti)/pz->mt->band->EffecHoleMass(Ti);
576 d_Fp_dpi += 2*e*mt->band->ThermalVp(Ti)*0.5*pcell->ilen[j]/pcell->area;
577 PetscScalar d_Fp_dpd = -2*e*pm*pz->mt->band->ThermalVp(Ti)*exp(-(Evi-Evj)/(k*Ti))*0.5*pcell->ilen[j]/pcell->area;
578 MatSetValue(*jac,zofs[z]+3*i+2,zofs[pz->pzone->zone_index]+3*n+2,-d_Fp_dpd,ADD_VALUES);
579 MatSetValue(*jac,N+zofs[z]+3*i+2,N+zofs[pz->pzone->zone_index]+3*n+2,-d_Fp_dpd,ADD_VALUES);
580 }
581 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
582 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
583 }
584 MatSetValues(*jac,3,IndexRe,3,ColRe,A.m,INSERT_VALUES);
585 MatSetValues(*jac,3,IndexIm,3,ColIm,A.m,INSERT_VALUES);
586 }
587
588 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A.m);
589 A=A/area;
590 if(pzone->zone_index < pz->pzone->zone_index)
591 {
592 A.m[0] = d_grad_P_dVi; //dfun(0)/dP(i)
593 A.m[1] = -e*pcell->area; //dfun(0)/dn(i)
594 A.m[2] = e*pcell->area; //dfun(0)/dp(i)
595 for(int j=0;j<ncell->nb_num;j++)
596 {
597 int nb = ncell->nb_array[j];
598 A.m[0] += -pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
599 PetscScalar value = pz->aux[n].eps*ncell->elen[j]/ncell->ilen[j];
600 MatSetValue(*jac,zofs[z]+3*i+0,zofs[pz->pzone->zone_index]+3*nb+0,value,INSERT_VALUES);
601 MatSetValue(*jac,N+zofs[z]+3*i+0,N+zofs[pz->pzone->zone_index]+3*nb+0,value,INSERT_VALUES);
602 }
603 MatSetValue(*jac,zofs[z]+3*i+0,zofs[pz->pzone->zone_index]+3*n+1,-e*ncell->area,INSERT_VALUES);
604 MatSetValue(*jac,zofs[z]+3*i+0,zofs[pz->pzone->zone_index]+3*n+2,e*ncell->area,INSERT_VALUES);
605 MatSetValue(*jac,N+zofs[z]+3*i+0,N+zofs[pz->pzone->zone_index]+3*n+1,-e*ncell->area,INSERT_VALUES);
606 MatSetValue(*jac,N+zofs[z]+3*i+0,N+zofs[pz->pzone->zone_index]+3*n+2,e*ncell->area,INSERT_VALUES);
607 }
608 else
609 {
610 A.m[0] = 1; //dfun(0)/dP(i)
611 A.m[1] = 0; //dfun(0)/dn(i)
612 A.m[2] = 0; //dfun(0)/dp(i)
613 MatSetValue(*jac,zofs[z]+3*i+0,zofs[pz->pzone->zone_index]+3*n+0,-1,INSERT_VALUES);
614 MatSetValue(*jac,N+zofs[z]+3*i+0,N+zofs[pz->pzone->zone_index]+3*n+0,-1,INSERT_VALUES);
615 }
616
617 A.m[4] += d_Fn_dni; //dfun(1)/dp(i)
618 A.m[8] += -d_Fp_dpi; //dfun(2)/dp(i)
619
620 MatSetValues(*jac,3,IndexRe,3,IndexRe,A.m,INSERT_VALUES);
621 MatSetValues(*jac,3,IndexIm,3,IndexIm,A.m,INSERT_VALUES);
622
623 Set_Mat3_zero(A);
624 A.m[4] = -omega;
625 A.m[8] = -omega;
626 MatSetValues(*jac,3,IndexIm,3,IndexRe,A.m,INSERT_VALUES);
627 A = -1.0*A;
628 MatSetValues(*jac,3,IndexRe,3,IndexIm,A.m,INSERT_VALUES);
629 Set_Vec3_zero(b);
630 VecSetValues(*r,3,IndexRe,b.v,INSERT_VALUES);
631 VecSetValues(*r,3,IndexIm,b.v,INSERT_VALUES);
632 }
633
634
635 void SMCZone::AC1_om_electrode(int i,PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,
636 vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
637 {
638
639 int equ_num = 3;
640 int size = pzone->davcell.size();
641 int z = zone_index;
642 int N;
643 MatGetSize(*jtmp,&N,PETSC_NULL);
644 int bc_index = electrode[i];
645 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
646 Vec3 J1,J2;
647 complex <PetscScalar> A[3];
648 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
649 PetscScalar R = pbc->R;
650 PetscScalar C = pbc->C;
651 PetscScalar L = pbc->L;
652 complex <PetscScalar> Z1(R,omega*L);
653 complex <PetscScalar> Y2(0,omega*C);
654 PetscInt EquRe = zofs[zone_index]+equ_num*size+i;
655 PetscInt EquIm = N + zofs[zone_index]+equ_num*size+i;
656 complex <PetscScalar> I(0.0,1.0);
657
658 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
659 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
660 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
661 {
662 int node=bc[bc_index].psegment->node_array[j];
663 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
664 complex <PetscScalar> dJdisp_dVi=0,dJdisp_dVr=0;
665 IndexRe[0] = zofs[z]+3*node+0;
666 IndexRe[1] = zofs[z]+3*node+1;
667 IndexRe[2] = zofs[z]+3*node+2;
668 IndexIm[0] = N + zofs[z]+3*node+0;
669 IndexIm[1] = N + zofs[z]+3*node+1;
670 IndexIm[2] = N + zofs[z]+3*node+2;
671 for(int k=0;k<pcell->nb_num;k++)
672 {
673 int nb = pcell->nb_array[k];
674 ColRe[0] = zofs[z]+3*nb+0;
675 ColRe[1] = zofs[z]+3*nb+1;
676 ColRe[2] = zofs[z]+3*nb+2;
677 ColIm[0] = N + zofs[z]+3*nb+0;
678 ColIm[1] = N + zofs[z]+3*nb+1;
679 ColIm[2] = N + zofs[z]+3*nb+2;
680
681 MatGetValues(*jtmp,1,&IndexRe[1],3,ColRe,J1.v);
682 MatGetValues(*jtmp,1,&IndexRe[2],3,ColRe,J2.v);
683 //for displacement current
684 dJdisp_dVi += omega*aux[node].eps*pcell->elen[k]/pcell->ilen[k]*I;
685 dJdisp_dVr = -omega*aux[node].eps*pcell->elen[k]/pcell->ilen[k]*I;
686 A[0]=(J1.v[0]-J2.v[0]+dJdisp_dVr)*DeviceDepth*Z1;
687 A[1]=(J1.v[1]-J2.v[1])*DeviceDepth*Z1;
688 A[2]=(J1.v[2]-J2.v[2])*DeviceDepth*Z1;
689 J1.v[0] = A[0].real();
690 J1.v[1] = A[1].real();
691 J1.v[2] = A[2].real();
692 J2.v[0] = A[0].imag();
693 J2.v[1] = A[1].imag();
694 J2.v[2] = A[2].imag();
695 MatSetValues(*jac,1,&EquRe,3,ColRe,J1.v,ADD_VALUES);
696 MatSetValues(*jac,1,&EquRe,3,ColIm,(-1.0*J2).v,ADD_VALUES);
697 MatSetValues(*jac,1,&EquIm,3,ColRe,J2.v,ADD_VALUES);
698 MatSetValues(*jac,1,&EquIm,3,ColIm,J1.v,ADD_VALUES);
699 }
700
701 MatGetValues(*jtmp,1,&IndexRe[1],3,IndexRe,J1.v);
702 MatGetValues(*jtmp,1,&IndexRe[2],3,IndexRe,J2.v);
703 A[0]=(J1.v[0]-J2.v[0]+dJdisp_dVi)*DeviceDepth*Z1;
704 A[1]=(J1.v[1]-J2.v[1])*DeviceDepth*Z1;
705 A[2]=(J1.v[2]-J2.v[2])*DeviceDepth*Z1;
706 J1.v[0] = A[0].real();
707 J1.v[1] = A[1].real();
708 J1.v[2] = A[2].real();
709 J2.v[0] = A[0].imag();
710 J2.v[1] = A[1].imag();
711 J2.v[2] = A[2].imag();
712 MatSetValues(*jac,1,&EquRe,3,IndexRe,J1.v,ADD_VALUES);
713 MatSetValues(*jac,1,&EquRe,3,IndexIm,(-1.0*J2).v,ADD_VALUES);
714 MatSetValues(*jac,1,&EquIm,3,IndexRe,J2.v,ADD_VALUES);
715 MatSetValues(*jac,1,&EquIm,3,IndexIm,J1.v,ADD_VALUES);
716 }
717 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
718 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
719
720 complex <PetscScalar> K = Z1*Y2 + PetscScalar(1.0);
721 MatSetValue(*jac,EquRe,EquRe,K.real(),INSERT_VALUES);
722 MatSetValue(*jac,EquRe,EquIm,-K.imag(),INSERT_VALUES);
723 MatSetValue(*jac,EquIm,EquRe,K.imag(),INSERT_VALUES);
724 MatSetValue(*jac,EquIm,EquIm,K.real(),INSERT_VALUES);
725
726 VecSetValue(*r,EquRe,pbc->Vac, INSERT_VALUES);
727 VecSetValue(*r,EquIm,0, INSERT_VALUES);
728
729 }
730
731
732 void SMCZone::AC1_stk_electrode(int i,PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,
733 vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
734 {
735 int equ_num = 3;
736 int size = pzone->davcell.size();
737 int z = zone_index;
738 int N;
739 MatGetSize(*jtmp,&N,PETSC_NULL);
740 int bc_index = electrode[i];
741 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
742 PetscScalar R = pbc->R;
743 PetscScalar C = pbc->C;
744 PetscScalar L = pbc->L;
745 complex <PetscScalar> Z1(R,omega*L);
746 complex <PetscScalar> Y2(0,omega*C);
747 complex <PetscScalar> I(0,1);
748 Vec3 J1,J2;
749 complex <PetscScalar> A[3];
750 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
751 PetscInt EquRe = zofs[zone_index]+equ_num*size+i;
752 PetscInt EquIm = N + zofs[zone_index]+equ_num*size+i;
753
754 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
755 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
756 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
757 {
758 int node=bc[bc_index].psegment->node_array[j];
759 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
760 PetscScalar VB = pbc->WorkFunction-aux[node].affinity;
761 PetscScalar ni = x[zofs[zone_index]+3*node+1]; //electron density of node i
762 PetscScalar pi = x[zofs[zone_index]+3*node+2]; //hole density of node i
763 PetscScalar dJ_dni = -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB)*0.5*pcell->ilen[0]
764 -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB)*0.5*pcell->ilen[pcell->nb_num-1];
765 PetscScalar dJ_dpi = -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB)*0.5*pcell->ilen[0]
766 -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB)*0.5*pcell->ilen[pcell->nb_num-1];
767 complex <PetscScalar> dJdisp_dVi=0,dJdisp_dVr=0;
768 IndexRe[0] = zofs[z]+3*node+0;
769 IndexRe[1] = zofs[z]+3*node+1;
770 IndexRe[2] = zofs[z]+3*node+2;
771 IndexIm[0] = N + zofs[z]+3*node+0;
772 IndexIm[1] = N + zofs[z]+3*node+1;
773 IndexIm[2] = N + zofs[z]+3*node+2;
774 for(int k=0;k<pcell->nb_num;k++)
775 {
776 int nb = pcell->nb_array[k];
777 ColRe[0] = zofs[z]+3*nb+0;
778 ColRe[1] = zofs[z]+3*nb+1;
779 ColRe[2] = zofs[z]+3*nb+2;
780 ColIm[0] = N + zofs[z]+3*nb+0;
781 ColIm[1] = N + zofs[z]+3*nb+1;
782 ColIm[2] = N + zofs[z]+3*nb+2;
783
784 //for displacement current
785 dJdisp_dVi += omega*aux[node].eps*pcell->elen[k]/pcell->ilen[k]*I;
786 dJdisp_dVr = -omega*aux[node].eps*pcell->elen[k]/pcell->ilen[k]*I;
787
788 A[0]=dJdisp_dVr*DeviceDepth*Z1;
789 A[1]=0*DeviceDepth*Z1;
790 A[2]=0*DeviceDepth*Z1;
791 J1.v[0] = A[0].real();
792 J1.v[1] = A[1].real();
793 J1.v[2] = A[2].real();
794 J2.v[0] = A[0].imag();
795 J2.v[1] = A[1].imag();
796 J2.v[2] = A[2].imag();
797 MatSetValues(*jac,1,&EquRe,3,ColRe,J1.v,ADD_VALUES);
798 MatSetValues(*jac,1,&EquRe,3,ColIm,(-1.0*J2).v,ADD_VALUES);
799 MatSetValues(*jac,1,&EquIm,3,ColRe,J2.v,ADD_VALUES);
800 MatSetValues(*jac,1,&EquIm,3,ColIm,J1.v,ADD_VALUES);
801 }
802 A[0] = dJdisp_dVi*DeviceDepth*Z1;
803 A[1] = dJ_dni*DeviceDepth*Z1;
804 A[2] = dJ_dpi*DeviceDepth*Z1;
805 J1.v[0] = A[0].real();
806 J1.v[1] = A[1].real();
807 J1.v[2] = A[2].real();
808 J2.v[0] = A[0].imag();
809 J2.v[1] = A[1].imag();
810 J2.v[2] = A[2].imag();
811 MatSetValues(*jac,1,&EquRe,3,IndexRe,J1.v,ADD_VALUES);
812 MatSetValues(*jac,1,&EquRe,3,IndexIm,(-1.0*J2).v,ADD_VALUES);
813 MatSetValues(*jac,1,&EquIm,3,IndexRe,J2.v,ADD_VALUES);
814 MatSetValues(*jac,1,&EquIm,3,IndexIm,J1.v,ADD_VALUES);
815 }
816 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
817 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
818
819 complex <PetscScalar> K = Z1*Y2 + PetscScalar(1.0);
820 MatSetValue(*jac,EquRe,EquRe,K.real(),INSERT_VALUES);
821 MatSetValue(*jac,EquRe,EquIm,-K.imag(),INSERT_VALUES);
822 MatSetValue(*jac,EquIm,EquRe,K.imag(),INSERT_VALUES);
823 MatSetValue(*jac,EquIm,EquIm,K.real(),INSERT_VALUES);
824
825 VecSetValue(*r,EquRe,pbc->Vac, INSERT_VALUES);
826 VecSetValue(*r,EquIm,0, INSERT_VALUES);
827 }
828
829
830 void SMCZone::AC1_ins_electrode(int i,PetscScalar omega, PetscScalar *x,Mat *jac,Vec *r,Mat *jtmp,
831 vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
832 {
833 int equ_num = 3;
834 int size = pzone->davcell.size();
835 int z = zone_index;
836 int N;
837 MatGetSize(*jtmp,&N,PETSC_NULL);
838 int bc_index = electrode[i];
839 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
840 PetscScalar R = pbc->R;
841 PetscScalar C = pbc->C;
842 PetscScalar L = pbc->L;
843
844 complex <PetscScalar> Z1(R,omega*L);
845 complex <PetscScalar> Y2(0,omega*C);
846 complex <PetscScalar> I(0,1);
847 Vec3 J1,J2;
848 complex <PetscScalar> A[3];
849 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
850 PetscInt EquRe = zofs[zone_index]+equ_num*size+i;
851 PetscInt EquIm = N + zofs[zone_index]+equ_num*size+i;
852
853 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
854 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
855 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
856 {
857 int node=bc[bc_index].psegment->node_array[j];
858 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
859 complex <PetscScalar> dJdisp_dVi=0,dJdisp_dVr=0;
860 IndexRe[0] = zofs[z]+3*node+0;
861 IndexRe[1] = zofs[z]+3*node+1;
862 IndexRe[2] = zofs[z]+3*node+2;
863 IndexIm[0] = N + zofs[z]+3*node+0;
864 IndexIm[1] = N + zofs[z]+3*node+1;
865 IndexIm[2] = N + zofs[z]+3*node+2;
866 for(int k=0;k<pcell->nb_num;k++)
867 {
868 int nb = pcell->nb_array[k];
869 ColRe[0] = zofs[z]+3*nb+0;
870 ColRe[1] = zofs[z]+3*nb+1;
871 ColRe[2] = zofs[z]+3*nb+2;
872 ColIm[0] = N + zofs[z]+3*nb+0;
873 ColIm[1] = N + zofs[z]+3*nb+1;
874 ColIm[2] = N + zofs[z]+3*nb+2;
875
876 //for displacement current
877 dJdisp_dVi += omega*aux[node].eps*pcell->elen[k]/pcell->ilen[k]*I;
878 dJdisp_dVr = -omega*aux[node].eps*pcell->elen[k]/pcell->ilen[k]*I;
879
880 A[0]=dJdisp_dVr*DeviceDepth*Z1;
881 A[1]=0.0;
882 A[2]=0.0;
883 J1.v[0] = A[0].real();
884 J1.v[1] = A[1].real();
885 J1.v[2] = A[2].real();
886 J2.v[0] = A[0].imag();
887 J2.v[1] = A[1].imag();
888 J2.v[2] = A[2].imag();
889 MatSetValues(*jac,1,&EquRe,3,ColRe,J1.v,ADD_VALUES);
890 MatSetValues(*jac,1,&EquRe,3,ColIm,(-1.0*J2).v,ADD_VALUES);
891 MatSetValues(*jac,1,&EquIm,3,ColRe,J2.v,ADD_VALUES);
892 MatSetValues(*jac,1,&EquIm,3,ColIm,J1.v,ADD_VALUES);
893 }
894 A[0] = dJdisp_dVi*DeviceDepth*Z1;
895 A[1] = 0.0;
896 A[2] = 0.0;
897 J1.v[0] = A[0].real();
898 J1.v[1] = A[1].real();
899 J1.v[2] = A[2].real();
900 J2.v[0] = A[0].imag();
901 J2.v[1] = A[1].imag();
902 J2.v[2] = A[2].imag();
903 MatSetValues(*jac,1,&EquRe,3,IndexRe,J1.v,ADD_VALUES);
904 MatSetValues(*jac,1,&EquRe,3,IndexIm,(-1.0*J2).v,ADD_VALUES);
905 MatSetValues(*jac,1,&EquIm,3,IndexRe,J2.v,ADD_VALUES);
906 MatSetValues(*jac,1,&EquIm,3,IndexIm,J1.v,ADD_VALUES);
907 }
908 MatAssemblyBegin(*jac,MAT_FLUSH_ASSEMBLY);
909 MatAssemblyEnd(*jac,MAT_FLUSH_ASSEMBLY);
910
911 complex <PetscScalar> K = Z1*Y2 + PetscScalar(1.0);
912 MatSetValue(*jac,EquRe,EquRe,K.real(),INSERT_VALUES);
913 MatSetValue(*jac,EquRe,EquIm,-K.imag(),INSERT_VALUES);
914 MatSetValue(*jac,EquIm,EquRe,K.imag(),INSERT_VALUES);
915 MatSetValue(*jac,EquIm,EquIm,K.real(),INSERT_VALUES);
916
917 VecSetValue(*r,EquRe,pbc->Vac, INSERT_VALUES);
918 VecSetValue(*r,EquIm,0, INSERT_VALUES);
919 }
920
921
922 void SMCZone::AC1_om_electrode_current(int i,PetscScalar omega, PetscScalar *x,Vec *v,Mat *jtmp,
923 vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth, PetscScalar &IacRe, PetscScalar &IacIm)
924 {
925 int equ_num = 3;
926 int size = pzone->davcell.size();
927 int z = zone_index;
928 int N;
929 MatGetSize(*jtmp,&N,PETSC_NULL);
930 int bc_index = electrode[i];
931 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(bc_index));
932 Mat3 A;
933 Vec3 V1,V2,b1,b2;
934 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
935 IacRe = 0;
936 IacIm = 0;
937 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
938 {
939 int node=bc[bc_index].psegment->node_array[j];
940 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
941
942 IndexRe[0] = zofs[z]+3*node+0;
943 IndexRe[1] = zofs[z]+3*node+1;
944 IndexRe[2] = zofs[z]+3*node+2;
945 IndexIm[0] = N + zofs[z]+3*node+0;
946 IndexIm[1] = N + zofs[z]+3*node+1;
947 IndexIm[2] = N + zofs[z]+3*node+2;
948 PetscScalar ViRe,ViIm;
949 VecGetValues(*v,1,&IndexRe[0],&ViRe);
950 VecGetValues(*v,1,&IndexIm[0],&ViIm);
951
952 for(int k=0;k<pcell->nb_num;k++)
953 {
954 int nb = pcell->nb_array[k];
955 ColRe[0] = zofs[z]+3*nb+0;
956 ColRe[1] = zofs[z]+3*nb+1;
957 ColRe[2] = zofs[z]+3*nb+2;
958 ColIm[0] = N + zofs[z]+3*nb+0;
959 ColIm[1] = N + zofs[z]+3*nb+1;
960 ColIm[2] = N + zofs[z]+3*nb+2;
961 MatGetValues(*jtmp,3,IndexRe,3,ColRe,A.m);
962 VecGetValues(*v,3,ColRe,V1.v);
963 VecGetValues(*v,3,ColIm,V2.v);
964 b1 = A*V1;
965 b2 = A*V2;
966 IacRe += (b1.v[1]-b1.v[2])*DeviceDepth;
967 IacIm += (b2.v[1]-b2.v[2])*DeviceDepth;
968 //for displacement current
969 PetscScalar VjRe = V1.v[0];
970 PetscScalar VjIm = V2.v[0];
971 IacRe += -omega*aux[node].eps*(ViIm-VjIm)/pcell->ilen[k]*pcell->elen[k]*DeviceDepth;
972 IacIm += omega*aux[node].eps*(ViRe-VjRe)/pcell->ilen[k]*pcell->elen[k]*DeviceDepth;
973 }
974 MatGetValues(*jtmp,3,IndexRe,3,IndexRe,A.m);
975 VecGetValues(*v,3,IndexRe,V1.v);
976 VecGetValues(*v,3,IndexIm,V2.v);
977 b1 = A*V1;
978 b2 = A*V2;
979 IacRe += (b1.v[1]-b1.v[2])*DeviceDepth;
980 IacIm += (b2.v[1]-b2.v[2])*DeviceDepth;
981 }
982 }
983
984
985 void SMCZone::AC1_stk_electrode_current(int i,PetscScalar omega, PetscScalar *x, Vec *v,Mat *jtmp,
986 vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth, PetscScalar &IacRe, PetscScalar &IacIm)
987 {
988 int equ_num = 3;
989 int size = pzone->davcell.size();
990 int z = zone_index;
991 int N;
992 MatGetSize(*jtmp,&N,PETSC_NULL);
993 int bc_index = electrode[i];
994 SchottkyBC *pbc = dynamic_cast <SchottkyBC * > (bc.Get_pointer(bc_index));
995 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
996 IacRe = 0;
997 IacIm = 0;
998 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
999 {
1000 int node=bc[bc_index].psegment->node_array[j];
1001 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1002 PetscScalar VB = pbc->WorkFunction-aux[node].affinity;
1003 PetscScalar ni = x[zofs[zone_index]+3*node+1]; //electron density of node i
1004 PetscScalar pi = x[zofs[zone_index]+3*node+2]; //hole density of node i
1005 PetscScalar dJ_dni = -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB)*0.5*pcell->ilen[0]
1006 -mt->band->pdSchottyJsn_pdn(ni,fs[node].T,VB)*0.5*pcell->ilen[pcell->nb_num-1];
1007 PetscScalar dJ_dpi = -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB)*0.5*pcell->ilen[0]
1008 -mt->band->pdSchottyJsp_pdp(pi,fs[node].T,VB)*0.5*pcell->ilen[pcell->nb_num-1];
1009
1010 IndexRe[0] = zofs[z]+3*node+0;
1011 IndexRe[1] = zofs[z]+3*node+1;
1012 IndexRe[2] = zofs[z]+3*node+2;
1013 IndexIm[0] = N + zofs[z]+3*node+0;
1014 IndexIm[1] = N + zofs[z]+3*node+1;
1015 IndexIm[2] = N + zofs[z]+3*node+2;
1016 PetscScalar ViRe,ViIm;
1017 PetscScalar niRe,niIm;
1018 PetscScalar piRe,piIm;
1019 VecGetValues(*v,1,&IndexRe[0],&ViRe);
1020 VecGetValues(*v,1,&IndexIm[0],&ViIm);
1021 VecGetValues(*v,1,&IndexRe[1],&niRe);
1022 VecGetValues(*v,1,&IndexIm[1],&niIm);
1023 VecGetValues(*v,1,&IndexRe[2],&piRe);
1024 VecGetValues(*v,1,&IndexIm[2],&piIm);
1025 for(int k=0;k<pcell->nb_num;k++)
1026 {
1027 int nb = pcell->nb_array[k];
1028 ColRe[0] = zofs[z]+3*nb+0;
1029 ColRe[1] = zofs[z]+3*nb+1;
1030 ColRe[2] = zofs[z]+3*nb+2;
1031 ColIm[0] = N + zofs[z]+3*nb+0;
1032 ColIm[1] = N + zofs[z]+3*nb+1;
1033 ColIm[2] = N + zofs[z]+3*nb+2;
1034
1035 //for displacement current
1036 PetscScalar VjRe,VjIm;
1037 VecGetValues(*v,1,&ColRe[0],&VjRe);
1038 VecGetValues(*v,1,&ColIm[0],&VjIm);
1039 IacRe += -omega*aux[node].eps*(ViIm-VjIm)/pcell->ilen[k]*pcell->elen[k]*DeviceDepth;
1040 IacIm += omega*aux[node].eps*(ViRe-VjRe)/pcell->ilen[k]*pcell->elen[k]*DeviceDepth;
1041 }
1042 IacRe += (dJ_dni*niRe + dJ_dpi*piRe)*DeviceDepth;
1043 IacIm += (dJ_dni*niIm + dJ_dpi*piIm)*DeviceDepth;
1044 }
1045 }
1046
1047 void SMCZone::AC1_ins_electrode_current(int i,PetscScalar omega, PetscScalar *x, Vec *v,Mat *jtmp,
1048 vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth, PetscScalar &IacRe, PetscScalar &IacIm)
1049 {
1050 int equ_num = 3;
1051 int size = pzone->davcell.size();
1052 int z = zone_index;
1053 int N;
1054 MatGetSize(*jtmp,&N,PETSC_NULL);
1055 int bc_index = electrode[i];
1056 InsulatorContactBC *pbc = dynamic_cast <InsulatorContactBC * > (bc.Get_pointer(bc_index));
1057 PetscInt IndexRe[3],IndexIm[3],ColRe[3],ColIm[3];
1058 IacRe = 0;
1059 IacIm = 0;
1060 //for displacement current
1061 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
1062 {
1063 int node=bc[bc_index].psegment->node_array[j];
1064 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
1065 IndexRe[0] = zofs[z]+3*node+0;
1066 IndexRe[1] = zofs[z]+3*node+1;
1067 IndexRe[2] = zofs[z]+3*node+2;
1068 IndexIm[0] = N + zofs[z]+3*node+0;
1069 IndexIm[1] = N + zofs[z]+3*node+1;
1070 IndexIm[2] = N + zofs[z]+3*node+2;
1071 PetscScalar ViRe,ViIm;
1072 VecGetValues(*v,1,&IndexRe[0],&ViRe);
1073 VecGetValues(*v,1,&IndexIm[0],&ViIm);
1074 for(int k=0;k<pcell->nb_num;k++)
1075 {
1076 int nb = pcell->nb_array[k];
1077 ColRe[0] = zofs[z]+3*nb+0;
1078 ColRe[1] = zofs[z]+3*nb+1;
1079 ColRe[2] = zofs[z]+3*nb+2;
1080 ColIm[0] = N + zofs[z]+3*nb+0;
1081 ColIm[1] = N + zofs[z]+3*nb+1;
1082 ColIm[2] = N + zofs[z]+3*nb+2;
1083
1084 PetscScalar VjRe,VjIm;
1085 VecGetValues(*v,1,&ColRe[0],&VjRe);
1086 VecGetValues(*v,1,&ColIm[0],&VjIm);
1087 IacRe += -omega*aux[node].eps*(ViIm-VjIm)/pcell->ilen[k]*pcell->elen[k]*DeviceDepth;
1088 IacIm += omega*aux[node].eps*(ViRe-VjRe)/pcell->ilen[k]*pcell->elen[k]*DeviceDepth;
1089 }
1090 }
1091 }
General-purpose Semiconductor Simulator