?  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 / ddm_nt1ac.cc
blobf4046b4c330eaf4f55aecbad8f996d2679d179b6
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 "mathfunc.h"
22 #include "ddm_nt1ac.h"
23 #include "log.h"
24
25 /* ----------------------------------------------------------------------------
26 * DDM_Solver_L1AC::form_linear_system_ac1: This function build the linear system.
27 */
28 void DDM_Solver_L1AC::form_linear_system_ac1(PetscScalar omega, PetscScalar *x, Mat *jac, Vec *r, Mat *jtmp)
29 {
30 MatZeroEntries(*jac);
31 VecZeroEntries(*r);
32
33 //Compute complex Jacobian matrix and right vector.
34 for(int z=0;z<zone_num;z++)
35 {
36 if(zonedata[z]->material_type==Semiconductor)
37 {
38 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
39
40 // process electrode.
41 for(int i=0;i<pzonedata->electrode.size();i++)
42 {
43 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
44 pzonedata->AC1_om_electrode(i,omega,x,jac,r,&JTmp,zofs,bc,DeviceDepth);
45 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
46 pzonedata->AC1_stk_electrode(i,omega,x,jac,r,&JTmp,zofs,bc,DeviceDepth);
47 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
48 pzonedata->AC1_ins_electrode(i,omega,x,jac,r,&JTmp,zofs,bc,DeviceDepth);
49 }
50
51 // process cell variables and boundaries.
52 for(int i=0;i<zone[z].davcell.size();i++)
53 {
54 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
55 SemiData *pfs = &pzonedata->fs[i];
56 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
57 pzonedata->AC1_ddm_inner(i,omega,x,jac,r,&JTmp,zofs);
58 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
59 pzonedata->AC1_ddm_ombc(i,omega,x,jac,r,&JTmp,zofs,bc);
60 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
61 pzonedata->AC1_ddm_stkbc(i,omega,x,jac,r,&JTmp,zofs,bc);
62 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
63 pzonedata->AC1_ddm_insulator_gate(i,omega,x,jac,r,&JTmp,zofs,bc);
64 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
65 {
66 InsulatorInterfaceBC *pbc;
67 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
68 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
69 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
70 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
71 pzonedata->AC1_ddm_interface(i,omega,x,jac,r,&JTmp,zofs,bc,pz,n_node);
72 }
73 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HomoInterface)
74 {
75 HomoInterfaceBC *pbc;
76 pbc = dynamic_cast<HomoInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
77 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
78 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
79 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
80 pzonedata->AC1_ddm_homojunction(i,omega,x,jac,r,&JTmp,zofs,bc,pz,n_node);
81 }
82 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HeteroInterface)
83 {
84 HeteroInterfaceBC *pbc;
85 pbc = dynamic_cast<HeteroInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
86 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
87 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
88 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
89 pzonedata->AC1_ddm_heterojunction(i,omega,x,jac,r,&JTmp,zofs,bc,pz,n_node);
90 }
91 else
92 pzonedata->AC1_ddm_inner(i,omega,x,jac,r,&JTmp,zofs);
93 }
94 }
95
96 // Insulator zone
97 if(zonedata[z]->material_type==Insulator)
98 {
99 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
100 for(int i=0;i<zone[z].davcell.size();i++)
101 {
102 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
103 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
104 pzonedata->AC1_ddm_inner(i,omega,x,r,jac,zofs);
105 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
106 pzonedata->AC1_ddm_gatebc(i,omega,x,r,jac,zofs,bc);
107 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==ChargedContact)
108 pzonedata->AC1_ddm_chargebc(i,omega,x,r,jac,zofs,bc);
109 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
110 {
111 ElectrodeInsulatorBC *pbc;
112 pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
113 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
114 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
115 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
116 pzonedata->AC1_ddm_electrode_insulator_interface(i,omega,x,r,jac,zofs,bc,pz,n_node);
117 }
118 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
119 {
120 InsulatorInterfaceBC *pbc;
121 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
122 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
123 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
124 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
125 pzonedata->AC1_ddm_semiconductor_insulator_interface(i,omega,x,r,jac,zofs,bc,pz,n_node);
126 }
127 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
128 {
129 InsulatorInsulatorBC *pbc;
130 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
131 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
132 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
133 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
134 pzonedata->AC1_ddm_insulator_insulator_interface(i,omega,x,r,jac,zofs,bc,pz,n_node);
135 }
136 else
137 pzonedata->AC1_ddm_inner(i,omega,x,r,jac,zofs);
138 }
139 for(int i=0;i<pzonedata->electrode.size();i++)
140 {
141 if(bc[pzonedata->electrode[i]].BCType==GateContact)
142 pzonedata->AC1_gate_electrode(i,omega,x,r,jac,zofs,bc,DeviceDepth);
143 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
144 pzonedata->AC1_charge_electrode(i,omega,x,r,jac,zofs,bc);
145 }
146 }
147
148 // Electrode zone
149 if(zonedata[z]->material_type==Conductor)
150 {
151 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
152 for(int i=0;i<zone[z].davcell.size();i++)
153 {
154 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
155 if(!pcell->bc_index || bc[pcell->bc_index-1].BCType==NeumannBoundary)
156 pzonedata->AC1_ddm_inner(i,omega,x,r,jac,zofs);
157 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
158 pzonedata->AC1_ddm_inner(i,omega,x,r,jac,zofs);
159 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==OhmicContact)
160 {
161 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
162 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
163 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
164 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
165 pzonedata->AC1_ddm_om_contact(i,omega,x,r,jac,zofs,bc,pz,n_node);
166 }
167 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==SchottkyContact)
168 {
169 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
170 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
171 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
172 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
173 pzonedata->AC1_ddm_stk_contact(i,omega,x,r,jac,zofs,bc,pz,n_node);
174 }
175 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==GateContact)
176 {
177 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
178 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
179 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
180 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
181 pzonedata->AC1_ddm_gate_contact(i,omega,x,r,jac,zofs,bc,pz,n_node);
182 }
183 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
184 {
185 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
186 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
187 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
188 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
189 pzonedata->AC1_ddm_charge_contact(i,omega,x,r,jac,zofs,bc,pz,n_node);
190 }
191 else
192 pzonedata->AC1_ddm_inner(i,omega,x,r,jac,zofs);
193 }
194 }
195
196 }
197 MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
198 MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
199 VecAssemblyBegin(*r);
200 VecAssemblyEnd(*r);
201 }
202
203
204 /* ----------------------------------------------------------------------------
205 * DDM_Solver_L1AC::Get_Jacobian: This function get a real Jacobian matrix.
206 */
207 void DDM_Solver_L1AC::Get_Jacobian()
208 {
209 // compute the scale of JTmp
210 zofs.resize(zone_num+1);
211 zofs[0] = 0;
212 N = 0;
213 for(int i=0;i<zone_num;i++)
214 {
215 if(zonedata[i]->material_type==Semiconductor) //semiconductor zone
216 {
217 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
218 N += 3*pzonedata->pzone->davcell.size();
219 N += pzonedata->electrode.size(); //additional equation for electrode
220 zofs[i+1] = N;
221 }
222 else if(zonedata[i]->material_type==Insulator) //Insulator zone
223 {
224 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
225 N += pzonedata->pzone->davcell.size();
226 N += pzonedata->electrode.size(); //additional equation for electrode
227 zofs[i+1] = N;
228 }
229 else if(zonedata[i]->material_type==Conductor) //Electrode zone
230 {
231 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
232 N += pzonedata->pzone->davcell.size();
233 zofs[i+1] = N;
234 }
235 else //other zones, such as PML and vacuum
236 {
237 zofs[i+1] = N;
238 }
239 }
240 // create DC solution array v
241 x = new PetscScalar[N];
242 // fill initial value
243 int offset=0;
244 for(int z=0;z<zone_num;z++)
245 {
246 if(zonedata[z]->material_type==Semiconductor)
247 {
248 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
249 for(int i=0;i<zone[z].davcell.size();i++)
250 {
251 x[offset+0] = pzonedata->fs[i].P;
252 x[offset+1] = pzonedata->fs[i].n;
253 x[offset+2] = pzonedata->fs[i].p;
254 offset += 3;
255 }
256 for(int j=0;j<pzonedata->electrode.size();j++)
257 {
258 PetscScalar P = bc.Get_pointer(pzonedata->electrode[j])->Get_Potential();
259 x[offset+0] = P;
260 offset += 1;
261 }
262 }
263
264 if(zonedata[z]->material_type==Insulator)
265 {
266 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
267 for(int i=0;i<zone[z].davcell.size();i++)
268 {
269 x[offset+0] = pzonedata->fs[i].P;
270 offset += 1;
271 }
272 for(int j=0;j<pzonedata->electrode.size();j++)
273 {
274 PetscScalar P = bc.Get_pointer(pzonedata->electrode[j])->Get_Potential();
275 x[offset+0] = P;
276 offset += 1;
277 }
278 }
279 if(zonedata[z]->material_type==Conductor)
280 {
281 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
282 for(int i=0;i<zone[z].davcell.size();i++)
283 {
284 x[offset+0] = pzonedata->fs[i].P;
285 offset += 1;
286 }
287 }
288 }
289 // Create Jacobian matrix data structure
290 // pre-alloc approximate memory
291 int *nnz = new int[N];
292 int *p = nnz;
293 for(int i=0;i<zone_num;i++)
294 {
295 //semiconductor zone
296 if(zonedata[i]->material_type==Semiconductor)
297 {
298 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
299 for(int j=0;j<pzonedata->pzone->davcell.size();j++)
300 {
301 const VoronoiCell* pcell = pzonedata->pzone->davcell.GetPointer(j);
302 //inner node, alloc exact memory.
303 if(!pcell->bc_index || bc[pcell->bc_index-1].psegment->interface==-1)
304 {
305 *p++ = 3*(pzonedata->pzone->davcell[j].nb_num+1);
306 *p++ = 3*(pzonedata->pzone->davcell[j].nb_num+1);
307 *p++ = 3*(pzonedata->pzone->davcell[j].nb_num+1);
308 }
309 else //interface node, slightly more than needed.
310 {
311 *p++ = 3*(2*pzonedata->pzone->davcell[j].nb_num+4);
312 *p++ = 3*(2*pzonedata->pzone->davcell[j].nb_num+4);
313 *p++ = 3*(2*pzonedata->pzone->davcell[j].nb_num+4);
314 }
315 }
316 for(int j=0;j<pzonedata->electrode.size();j++)
317 {
318 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
319 }
320
321 }
322 //Insulator zones, poisson equation only
323 if(zonedata[i]->material_type==Insulator)
324 {
325 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
326 for(int j=0;j<zone[i].davcell.size();j++)
327 {
328 *p++ = zone[i].davcell[j].nb_num+1;
329 }
330 for(int j=0;j<pzonedata->electrode.size();j++)
331 {
332 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
333 }
334 }
335 //Electrode zones, poisson equation only
336 if(zonedata[i]->material_type==Conductor)
337 {
338 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
339 for(int j=0;j<zone[i].davcell.size();j++)
340 *p++ = 1*(zone[i].davcell[j].nb_num+1);
341 }
342 }
343 MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,0,nnz,&JTmp);
344
345 //Compute Jacobian entries for flux along triangle edges.
346 for(int z=0;z<zone_num;z++)
347 {
348 if(zonedata[z]->material_type==Semiconductor)
349 {
350 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
351 // compute flux Jacobian along triangle edges
352 for(int i=0;i<zone[z].datri.size();i++)
353 {
354 Tri *ptri = &zone[z].datri[i];
355 pzonedata->J1E_Tri_ddm(ptri,x,&JTmp,zofs);
356 }
357 }
358 }
359 MatAssemblyBegin(JTmp,MAT_FINAL_ASSEMBLY);
360 MatAssemblyEnd(JTmp,MAT_FINAL_ASSEMBLY);
361 delete [] nnz;
362
363 }
364
365 /* ----------------------------------------------------------------------------
366 * DDM_Solver_L1AC::init_solver: This function do initial setup for AC linear solver
367 */
368 int DDM_Solver_L1AC::init_solver(SolveDefine &sv)
369 {
370 gss_log.string_buf()<<"DDM solver Level 1 AC Sweep init...";
371
372 //distable advaced features.
373 for(int i=0;i<zone_num;i++)
374 if(zonedata[i]->material_type == Semiconductor)
375 {
376 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
377 pzonedata->HighFieldMobility= sv.HighFieldMobility;
378 pzonedata->ImpactIonization = sv.ImpactIonization;
379 pzonedata->IIType = sv.IIType;
380 pzonedata->BandBandTunneling = sv.BandBandTunneling;
381 pzonedata->IncompleteIonization = sv.IncompleteIonization;
382 pzonedata->QuantumMechanical = sv.QuantumMechanical;
383 pzonedata->Fermi = sv.Fermi;
384 pzonedata->EJModel = sv.EJModel;
385 }
386
387 Get_Jacobian();
388
389 // scale of complex Jac
390 VecCreateSeq(PETSC_COMM_SELF,2*N,&v);
391 VecDuplicate(v,&r);
392
393 // Create Jacobian matrix data structure
394 // pre-alloc approximate memory
395 int *nnz = new int[2*N];
396 int *p = nnz;
397 //fill half of the array nnz
398 for(int i=0;i<zone_num;i++)
399 {
400 //semiconductor zone
401 if(zonedata[i]->material_type==Semiconductor)
402 {
403 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
404 for(int j=0;j<pzonedata->pzone->davcell.size();j++)
405 {
406 const VoronoiCell* pcell = pzonedata->pzone->davcell.GetPointer(j);
407 //inner node, alloc exact memory.
408 if(!pcell->bc_index || bc[pcell->bc_index-1].psegment->interface==-1)
409 {
410 *p++ = 2*3*(pzonedata->pzone->davcell[j].nb_num+1);
411 *p++ = 2*3*(pzonedata->pzone->davcell[j].nb_num+1);
412 *p++ = 2*3*(pzonedata->pzone->davcell[j].nb_num+1);
413 }
414 else //interface node, slightly more than needed.
415 {
416 *p++ = 2*3*(2*pzonedata->pzone->davcell[j].nb_num+4);
417 *p++ = 2*3*(2*pzonedata->pzone->davcell[j].nb_num+4);
418 *p++ = 2*3*(2*pzonedata->pzone->davcell[j].nb_num+4);
419 }
420 }
421 for(int j=0;j<pzonedata->electrode.size();j++)
422 {
423 *p++ = 2*(bc[pzonedata->electrode[j]].psegment->node_array.size()+1);
424 }
425
426 }
427 //Insulator zones, poisson equation only
428 if(zonedata[i]->material_type==Insulator)
429 {
430 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
431 for(int j=0;j<zone[i].davcell.size();j++)
432 {
433 *p++ = 2*(zone[i].davcell[j].nb_num+1);
434 }
435 for(int j=0;j<pzonedata->electrode.size();j++)
436 {
437 *p++ = 2*(bc[pzonedata->electrode[j]].psegment->node_array.size()+1);
438 }
439 }
440 //Electrode zones, poisson equation only
441 if(zonedata[i]->material_type==Conductor)
442 {
443 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
444 for(int j=0;j<zone[i].davcell.size();j++)
445 *p++ = 2*(zone[i].davcell[j].nb_num+1);
446 }
447
448 }
449 //fill another half of nnz array
450 for(int i=0;i<N;i++)
451 nnz[N+i]=nnz[i];
452
453 MatCreate(PETSC_COMM_SELF,&J);
454 MatSetSizes(J,2*N,2*N,2*N,2*N);
455 if(sv.LS=="superlu")
456 {
457 #ifdef PETSC_HAVE_SUPERLU
458 MatSetType(J,MATSUPERLU);
459 #else
460 MatSetType(J,MATSEQAIJ);
461 #endif
462 }
463 else if(sv.LS=="umfpack")
464 {
465 #ifdef PETSC_HAVE_UMFPACK
466 MatSetType(J,MATUMFPACK);
467 #else
468 MatSetType(J,MATSEQAIJ);
469 #endif
470 }
471 else
472 {
473 MatSetType(J,MATSEQAIJ);
474 }
475 MatSeqAIJSetPreallocation(J,0,nnz);
476
477 //MatCreateSeqAIJ(PETSC_COMM_SELF,2*N,2*N,0,nnz,&J);
478 delete [] nnz;
479
480 KSPCreate(PETSC_COMM_SELF,&ksp);
481 KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN);
482 KSPGetPC(ksp,&pc);
483 //KSPSetMonitor(ksp,KSPDefaultMonitor,PETSC_NULL,PETSC_NULL);
484 if(sv.LS=="lu"||sv.LS=="superlu"||sv.LS=="umfpack")
485 {
486 KSPSetType(ksp,KSPGMRES);
487 KSPSetTolerances(ksp,1e-20*N,1e-18*N,PETSC_DEFAULT,N/2);
488 PCSetType(pc,PCLU);
489 PCFactorSetReuseOrdering(pc,PETSC_TRUE);
490 PCFactorSetPivoting(pc,1.0);
491 PCFactorReorderForNonzeroDiagonal(pc,1e-14);
492 PCFactorSetShiftNonzero(pc,1e-12);
493 }
494 else
495 {
496 KSPSetType(ksp,sv.LS.c_str());
497 if(sv.LS=="gmres") KSPGMRESSetRestart(ksp,150);
498 KSPSetTolerances(ksp,1e-20*N,1e-18*N,PETSC_DEFAULT,max(35,N/2));
499 PCSetType(pc,PCILU);
500 PCFactorSetLevels(pc,3);
501 PCFactorSetShiftNonzero(pc,1e-14);
502 PCFactorReorderForNonzeroDiagonal(pc,1e-12);
503 }
504 KSPSetFromOptions(ksp);
505 gss_log.string_buf()<<"done\n";
506 gss_log.record();
507 return 0;
508 }
509
510
511 /* ----------------------------------------------------------------------------
512 * DDM_Solver_L1AC::do_solve: This function solve the problem
513 */
514 int DDM_Solver_L1AC::do_solve(SolveDefine &sv)
515 {
516 if(sv.Type==ACSWEEP)
517 solve_ac(sv);
518 else
519 {
520 gss_log.string_buf()<<"DDML1AC not support this solver type!\n";
521 gss_log.record();
522 }
523 return 0;
524 }
525
526
527 /* ----------------------------------------------------------------------------
528 * DDM_Solver_L1AC::solve_ac: This function compute device response to
529 * small AC signal AC.
530 */
531 int DDM_Solver_L1AC::solve_ac(SolveDefine &sv)
532 {
533 FILE *fiv;
534 PetscScalar current_scale_mA = scale_unit.s_mA;
535 PetscScalar voltage_scale_V =scale_unit.s_volt;
536 KSPConvergedReason reason;
537 PetscInt its;
538 PetscReal rnorm;
539 bc.Update_Vapp(0);
540 bc.Update_Iapp(0);
541
542 //set which electrodes are required to record IV
543 for(int i=0;i<zone_num;i++)
544 {
545 if(zonedata[i]->material_type == Semiconductor)
546 {
547 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
548 pzonedata->ImpactIonization = sv.ImpactIonization;
549 pzonedata->IncompleteIonization = sv.IncompleteIonization;
550 pzonedata->QuantumMechanical = sv.QuantumMechanical;
551 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
552 for(int j=0;j<pzonedata->electrode.size();j++)
553 {
554 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
555 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
556 }
557 }
558 if(zonedata[i]->material_type == Insulator)
559 {
560 ISZone * pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
561 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
562 for(int j=0;j<pzonedata->electrode.size();j++)
563 {
564 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
565 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
566 }
567 }
568 }
569
570 // output DC Scan information
571 gss_log.string_buf()<<"AC scan from "<<sv.FStart*scale_unit.s_second/1e6
572 <<" MHz step size "<<sv.FMultiple
573 <<" to "<<sv.FStop*scale_unit.s_second/1e6<<" MHz"<<"\n";
574 gss_log.record();
575
576 // prepare IV file. If no file given, output to screen.
577 if(!sv.IVFile.empty())
578 {
579 fiv=fopen(sv.IVFile.c_str(),"w");
580 fprintf(fiv,"#");
581 fprintf(fiv,"Freq[MHz] ");
582 for(int j=0;j<sv.Electrode_Record.size();j++)
583 fprintf(fiv,"Vp(%s)[V] I(%s)[mA] ",
584 sv.Electrode_Record_Name[j].c_str(),
585 sv.Electrode_Record_Name[j].c_str()
586 );
587 fprintf(fiv,"\n");
588 }
589 else
590 fiv = stdout;
591
592 // begin
593 PetscScalar omega = 2*PI*sv.FStart;
594 if(omega<=0)
595 {
596 gss_log.string_buf()<<"The start frequency shouldn't be zero or negtive.\n";
597 gss_log.record();
598 return 0;
599 }
600
601 bc.Set_electrode_type(sv.Electrode_ACScan_Name.c_str(),VoltageBC);
602 bc.Set_Vac(sv.Electrode_ACScan_Name.c_str(),sv.VAC);
603
604 while( omega < 2*PI*sv.FStop)
605 {
606 gss_log.string_buf()<<"AC Scan: f("<<sv.Electrode_ACScan_Name<<") = "
607 <<omega/(2*PI)*scale_unit.s_second/1e6<<" MHz "<<"\n";
608 gss_log.record();
609 form_linear_system_ac1(omega,x,&J,&r,&JTmp);
610 KSPSetUp(ksp);
611 KSPSolve(ksp,r,v);
612 KSPGetConvergedReason(ksp,&reason);
613 KSPGetIterationNumber(ksp,&its);
614 KSPGetResidualNorm(ksp,&rnorm);
615 gss_log.string_buf()<<"------> residual norm = "<<rnorm<<" its = "<<its<<" with "<<SNESConvergedReasons[reason]<<"\n\n";
616 gss_log.record();
617 compute_electrode_current(omega);
618 fprintf(fiv,"%e\t",double(omega/(2*PI)*scale_unit.s_second/1e6));
619 for(int j=0;j<sv.Electrode_Record.size();j++)
620 {
621 int zone_index = bc[sv.Electrode_Record[j]].psegment->zone_index;
622 complex<PetscScalar> I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Iac()/current_scale_mA;
623 complex<PetscScalar> P = bc.Get_pointer(sv.Electrode_Record[j])->Get_Pac()/voltage_scale_V;
624 fprintf(fiv,"%e %e %e\t %e %e %e\t",
625 double(P.real()),double(P.imag()),double(abs(P)),
626 double(I.real()),double(I.imag()),double(abs(I)));
627 }
628 fprintf(fiv,"\n");
629 fflush(fiv);
630 omega*=sv.FMultiple;
631 }
632
633 if(!sv.IVFile.empty()) fclose(fiv);
634 return 0;
635 }
636
637
638 /* ----------------------------------------------------------------------------
639 * DDM_Solver_L1AC::compute_electrode_current: This function compute electrode current after
640 * AC sweep.
641 */
642 void DDM_Solver_L1AC::compute_electrode_current(PetscScalar omega)
643 {
644 //------------------------------------------------------------
645 PetscScalar *vv;
646 VecGetArray(v,&vv);
647
648 for(int z=0;z<zone_num;z++)
649 {
650 if(zonedata[z]->material_type==Semiconductor)
651 {
652 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
653 int equ_num = 3;
654 int size = pzonedata->pzone->davcell.size();
655 for(int i=0;i<pzonedata->electrode.size();i++)
656 {
657 PetscScalar IacRe=0,IacIm=0;
658 PetscScalar PRe = vv[zofs[z]+size*equ_num+i];
659 PetscScalar PIm = vv[N+zofs[z]+size*equ_num+i];
660
661 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
662 pzonedata->AC1_om_electrode_current(i,omega,x,&v,&JTmp, zofs, bc, DeviceDepth, IacRe, IacIm);
663 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
664 pzonedata->AC1_stk_electrode_current(i,omega,x,&v,&JTmp, zofs, bc, DeviceDepth, IacRe, IacIm);
665 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
666 pzonedata->AC1_ins_electrode_current(i,omega,x,&v,&JTmp, zofs, bc, DeviceDepth, IacRe, IacIm);
667 complex<PetscScalar> I(IacRe,IacIm);
668 complex<PetscScalar> P(PRe,PIm);
669 bc.Get_pointer(pzonedata->electrode[i])->Set_Iac(I);
670 bc.Get_pointer(pzonedata->electrode[i])->Set_Pac(P);
671
672 }
673 }
674 if(zonedata[z]->material_type==Insulator)
675 {
676 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
677 int equ_num = 1;
678 int size = pzonedata->pzone->davcell.size();
679 for(int i=0;i<pzonedata->electrode.size();i++)
680 {
681 PetscScalar IacRe=0,IacIm=0;
682 PetscScalar PRe = vv[zofs[z]+size*equ_num+i];
683 PetscScalar PIm = vv[N+zofs[z]+size*equ_num+i];
684 if(bc[pzonedata->electrode[i]].BCType==GateContact)
685 pzonedata->AC1_gate_electrode_current(i,omega,x,&v,zofs, bc, DeviceDepth, IacRe, IacIm);
686 complex<PetscScalar> I(IacRe,IacIm);
687 complex<PetscScalar> P(PRe,PIm);
688 bc.Get_pointer(pzonedata->electrode[i])->Set_Iac(I);
689 bc.Get_pointer(pzonedata->electrode[i])->Set_Pac(P);
690 }
691 }
692 }
693
694 VecRestoreArray(v,&vv);
695 }
696
697
698 /* ----------------------------------------------------------------------------
699 * DDM_Solver_L1AC::destroy_solver: This function do destroy the nonlinear solver
700 */
701 int DDM_Solver_L1AC::destroy_solver(SolveDefine &sv)
702 {
703 // free work space
704 N = 0;
705 zofs.clear();
706 delete [] x;
707 VecDestroy(v);
708 VecDestroy(r);
709 MatDestroy(J);
710 MatDestroy(JTmp);
711 KSPDestroy(ksp);
712 return 0;
713 }
714
General-purpose Semiconductor Simulator