?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
initial checkin, based on GSS 0.46 CVS
[gss-tcad.git] / src / solver / data_init.cc
blobb522a56afc171c9ecbaba3d74534b97e5f28c675
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: March 13, 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 "bsolver.h"
22 #include "log.h"
23 #include "typedef.h"
24 #include "adolc.h"
25 unsigned int adtl::AutoDScalar::numdir = 9;
26
27 int BSolver::build_zonedata()
28 {
29 zonedata.resize(zone_num);
30 // determine zone material by zone label. the label is compatible with medici
31 for(int i=0; i<zone_num; i++)
32 { //it is an Insulator zone
33 if(IsInsulator(zone[i].zonelabel))
34 {
35 zonedata[i] = new ISZone;
36 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
37 pzonedata->pzone = &zone[i];
38 pzonedata->zone_index = zone[i].zone_index;
39 pzonedata->pzone = &zone[i];
40 pzonedata->node_num = zone[i].davcell.size();
41 pzonedata->tri_num = zone[i].datri.size();
42 pzonedata->material_type = Insulator;
43 pzonedata->material = IsInsulator(zone[i].zonelabel);
44 try {pzonedata->mt = new MatInsulator(zone[i].zonename,
45 FormatInsulatorString(zone[i].zonelabel),&scale_unit);}
46 catch(int){return 1;}
47 pzonedata->fs = new ISData[pzonedata->node_num];
48 pzonedata->aux = new ISAuxData[pzonedata->node_num];
49 pzonedata->electrode.clear();
50 for(int j=0;j<pzonedata->pzone->dasegment.size();j++)
51 {
52 int bc_index = pzonedata->pzone->dasegment[j].bc_index-1;
53 if(bc[bc_index].BCType==GateContact ||
54 bc[bc_index].BCType==ChargedContact )
55 pzonedata->electrode.push_back(bc_index);
56 }
57 //pzonedata->report();
58 }
59 //if this is a semiconductor zone
60 else if(IsSemiconductor(zone[i].zonelabel))
61 {
62 zonedata[i] = new SMCZone;
63 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
64 pzonedata->pzone = &zone[i];
65 pzonedata->zone_index = zone[i].zone_index;
66 pzonedata->material = IsSemiconductor(zone[i].zonelabel);
67 pzonedata->material_type = Semiconductor;
68 try{pzonedata->mt = new MatSemiconductor(zone[i].zonename,
69 FormatSemiconductorString(zone[i].zonelabel),PMIS_define,&scale_unit);}
70 catch(int){return 1;}
71 pzonedata->node_num = zone[i].davcell.size();
72 pzonedata->tri_num = zone[i].datri.size();
73 pzonedata->fs = new SemiData[pzonedata->node_num];
74 pzonedata->aux = new SemiAuxData[pzonedata->node_num];
75 pzonedata->electrode.clear();
76 for(int j=0;j<pzonedata->pzone->dasegment.size();j++)
77 {
78 int bc_index = pzonedata->pzone->dasegment[j].bc_index-1;
79 if(bc[bc_index].BCType==OhmicContact||
80 bc[bc_index].BCType==SchottkyContact||
81 bc[bc_index].BCType==InsulatorContact)
82 pzonedata->electrode.push_back(bc_index);
83 }
84 //pzonedata->report();
85 }
86 //oh it is an electrode,medici can generate this useless region
87 else if(IsElectrode(zone[i].zonelabel))
88 {
89 zonedata[i] = new ElZone;
90 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
91 pzonedata->pzone = &zone[i];
92 pzonedata->zone_index = zone[i].zone_index;
93 pzonedata->node_num = zone[i].davcell.size();
94 pzonedata->tri_num = zone[i].datri.size();
95 pzonedata->material_type = Conductor;
96 pzonedata->material = Conductor;
97 try {pzonedata->mt = new MatConductor(zone[i].zonename,
98 FormatElectrodeString(zone[i].zonelabel),&scale_unit);}
99 catch(int){return 1;}
100 pzonedata->fs = new ELData[pzonedata->node_num];
101 pzonedata->aux = new ELAuxData[pzonedata->node_num];
102 pzonedata->electrode.clear();
103 for(int j=0;j<pzonedata->pzone->dasegment.size();j++)
104 {
105 int bc_index = pzonedata->pzone->dasegment[j].bc_index-1;
106 if(bc[bc_index].BCType==OhmicContact||
107 bc[bc_index].BCType==SchottkyContact||
108 bc[bc_index].BCType==GateContact)
109 pzonedata->electrode.push_back(bc_index);
110 }
111 if(pzonedata->electrode.size()>1)
112 {
113 gss_log.string_buf()<<"Electrode region "<<zone[i].zonename<<" has more than one Contact, not supported yet.\n";
114 gss_log.record();
115 return 1;
116 }
117 //pzonedata->report();
118 }
119 //Data for Vacuum zone
120 else if(IsVacuum(zone[i].zonelabel))
121 {
122 zonedata[i] = new VacuumZone;
123 VacuumZone *pzonedata = dynamic_cast< VacuumZone * >(zonedata[i]);
124 pzonedata->pzone = &zone[i];
125 pzonedata->zone_index = zone[i].zone_index;
126 pzonedata->node_num = zone[i].davcell.size();
127 pzonedata->tri_num = zone[i].datri.size();
128 pzonedata->material_type = Vacuum;
129 pzonedata->material = Vacuum;
130 try {pzonedata->mt = new MatVacuum(zone[i].zonename,
131 FormatVacuumString(zone[i].zonelabel),&scale_unit);}
132 catch(int){return 1;}
133 pzonedata->fs = new VacuumData[pzonedata->node_num];
134 pzonedata->aux = new VacuumAuxData[pzonedata->node_num];
135 }
136 //PML data
137 else if(IsPML(zone[i].zonelabel))
138 {
139 zonedata[i] = new PMLZone;
140 PMLZone *pzonedata = dynamic_cast< PMLZone * >(zonedata[i]);
141 pzonedata->pzone = &zone[i];
142 pzonedata->zone_index = zone[i].zone_index;
143 pzonedata->node_num = zone[i].davcell.size();
144 pzonedata->tri_num = zone[i].datri.size();
145 pzonedata->material_type = PML;
146 pzonedata->material = PML;
147 try {pzonedata->mt = new MatPML(zone[i].zonename,
148 FormatPMLString(zone[i].zonelabel),&scale_unit);}
149 catch(int){return 1;}
150 pzonedata->fs = new PMLData[pzonedata->node_num];
151 pzonedata->aux = new PMLAuxData[pzonedata->node_num];
152 }
153 else
154 {
155 gss_log.string_buf()<<"error, no such region type "<<zone[i].zonelabel<<" specified!\n";
156 gss_log.record();
157 return 1;
158 }
159 }
160
161 // process inter_connect here
162 for(int i=0; i<bc.size(); i++)
163 {
164 if(bc.Get_pointer(i)->BCType==OhmicContact)
165 {
166 OhmicBC *pbc = dynamic_cast <OhmicBC * > (bc.Get_pointer(i));
167 if(pbc->connect_elec!="")
168 {
169 int flag=0;
170 for(int j=0;j<bc.size();j++)
171 {
172 if(pbc->connect_elec==bc.Get_pointer(j)->psegment->label&&
173 IsSemiconductor(zone[bc.Get_pointer(j)->psegment->zone_index].zonelabel))
174 {
175 flag=1;
176 pbc->inner_connect=j;
177 pbc->connect_zone = bc.Get_pointer(j)->psegment->zone_index;
178 if(bc.Get_pointer(pbc->inner_connect)->BCType==OhmicContact)
179 {
180 OhmicBC *om_bc = dynamic_cast <OhmicBC * >(bc.Get_pointer(pbc->inner_connect));
181 om_bc->inner_connect=i;
182 om_bc->connect_zone=pbc->psegment->zone_index;
183 om_bc->connect_elec=pbc->psegment->label;
184 om_bc->pzonedata=zonedata[pbc->psegment->zone_index];
185 }
186 else
187 {
188 sprintf(log_buf,"\nCheck Inner Connect failed, BC %s should connect to an Ohmic/Schottky Electrode!\n",
189 bc.Get_pointer(pbc->inner_connect)->psegment->label);
190 GSS_LOG();
191 return 1;
192 }
193 }
194 }
195 if(!flag)
196 {
197 sprintf(log_buf,"\nSorry, I can't find a suitable inner conect for %s!\n",pbc->connect_elec.c_str());
198 GSS_LOG();
199 return 1;
200 }
201 }
202 }
203 }
204
205 return 0;
206 }
207
208 int BSolver::setup_init_data()
209 {
210 for(int i=0; i<zone_num; i++)
211 if(zonedata[i]->Init(&zone[i],LatticeTemp,&scale_unit)) return 1;
212 return 0;
213 }
214
215
216 int BSolver::setup_doping()
217 {
218 for(int z=0; z<zone_num; z++)
219 if(zonedata[z]->material_type == Semiconductor)
220 {
221 SMCZone *pzonedata = dynamic_cast<SMCZone *> (zonedata[z]);
222 for(int i=0;i<pzonedata->node_num;i++)
223 {
224 double x = pzonedata->pzone->danode[i].x;
225 double y = pzonedata->pzone->danode[i].y;
226 for(int j=0;j<doping_func.size();j++)
227 {
228 double doping = doping_func[j]->profile(x,y);
229 if(doping > 0)
230 pzonedata->aux[i].Nd += doping;
231 else
232 pzonedata->aux[i].Na += fabs(doping);
233 }
234 }
235 }
236 return 0;
237 }
238
239
240 int BSolver::import_doping_from_cgns(char* cgnsfile)
241 {
242 for(int z=0; z<zone_num; z++)
243 {
244 if(zonedata[z]->material_type == Semiconductor)
245 {
246 SMCZone *pzonedata = dynamic_cast<SMCZone *> (zonedata[z]);
247 if(pzonedata->import_doping(cgnsfile,&scale_unit)) return 1;
248 }
249 }
250 return 0;
251 }
252
253
254 int BSolver::import_mole_from_cgns(char* cgnsfile)
255 {
256 for(int z=0; z<zone_num; z++)
257 {
258 if(IsSingleCompSemiconductor(zone[z].zonelabel))
259 {
260 SMCZone *pzonedata = dynamic_cast<SMCZone *> (zonedata[z]);
261 if(pzonedata->import_mole(cgnsfile,1)) return 1;
262 }
263 }
264 return 0;
265 }
266
267 int BSolver::setup_bc()
268 {
269 //the bc init here
270 Segment **segment = new Segment*[bc_num];
271 for(int i=1; i<=bc_num;i++)
272 segment[i-1] = Get_segment(i);
273 if(bc.InitBC(bc_num,segment,zone_num,zone,zinterface,pcmdbuf,LatticeTemp,&scale_unit)) return 1;
274 delete [] segment;
275
276 //for bc marks are assigned to edge, some nodes may have two bc marks
277 //assign this type of nodes to correct bc type
278 for(int i=0;i<bc.size();i++)
279 {
280 if(bc[i].psegment->interface==-1)
281 {
282 int zone_index = bc[i].psegment->zone_index;
283 for(int j=0;j<bc[i].psegment->node_array.size();j++)
284 {
285 int node_index = bc[i].psegment->node_array[j];
286 int bc_index=zone[zone_index].danode[node_index].bc_index;
287 if(bc[bc_index-1].BCType < bc[i].BCType)
288 {
289 zone[zone_index].danode[node_index].bc_index = i+1;
290 zone[zone_index].davcell.GetPointer(node_index)->bc_index = i+1;
291 }
292 }
293 }
294 else
295 {
296 int inf = bc[i].psegment->interface;
297 int z1 = zinterface[inf].zone1;
298 int z2 = zinterface[inf].zone2;
299 if(z1==bc[i].psegment->zone_index)
300 for(int j=0;j<zinterface[inf].index_array1.size();j++)
301 {
302 int node_index = zinterface[inf].index_array1[j];
303 int bc_index=zone[z1].danode[node_index].bc_index;
304 if(bc[bc_index-1].BCType < bc[i].BCType)
305 {
306 zone[z1].danode[node_index].bc_index = i+1;
307 zone[z1].davcell.GetPointer(node_index)->bc_index = i+1;
308 }
309 }
310 if(z2==bc[i].psegment->zone_index)
311 for(int j=0;j<zinterface[inf].index_array2.size();j++)
312 {
313 int node_index = zinterface[inf].index_array2[j];
314 int bc_index=zone[z2].danode[node_index].bc_index;
315 if(bc[bc_index-1].BCType < bc[i].BCType)
316 {
317 zone[z2].danode[node_index].bc_index = i+1;
318 zone[z2].davcell.GetPointer(node_index)->bc_index = i+1;
319 }
320 }
321 }
322 }
323
324 return 0;
325 }
326
327
328 void BSolver::reorder()
329 {
330
331 for(int z=0; z<zone_num; z++)
332 {
333 int * old_to_new = reorder_zone(z);
334 if(zonedata[z]->material_type == Semiconductor)
335 {
336 SMCZone *pzonedata = dynamic_cast<SMCZone *> (zonedata[z]);
337 SemiData* fs_new = new SemiData[pzonedata->node_num];
338 SemiAuxData* aux_new = new SemiAuxData[pzonedata->node_num];
339 for(int i=0;i<pzonedata->node_num;i++)
340 {
341 fs_new[old_to_new[i]] = pzonedata->fs[i];
342 aux_new[old_to_new[i]] = pzonedata->aux[i];
343 }
344 delete [] pzonedata->fs;
345 delete [] pzonedata->aux;
346 pzonedata->fs = fs_new;
347 pzonedata->aux = aux_new;
348 }
349 if(zonedata[z]->material_type == Insulator)
350 {
351 ISZone *pzonedata = dynamic_cast<ISZone *> (zonedata[z]);
352 ISData* fs_new = new ISData[zonedata[z]->node_num];
353 ISAuxData* aux_new = new ISAuxData[pzonedata->node_num];
354 for(int i=0;i<pzonedata->node_num;i++)
355 {
356 fs_new[old_to_new[i]] = pzonedata->fs[i];
357 aux_new[old_to_new[i]] = pzonedata->aux[i];
358 }
359 delete [] pzonedata->fs;
360 delete [] pzonedata->aux;
361 pzonedata->fs = fs_new;
362 pzonedata->aux = aux_new;
363 }
364 if(zonedata[z]->material_type == Conductor)
365 {
366 ElZone *pzonedata = dynamic_cast<ElZone *> (zonedata[z]);
367 ELData* fs_new = new ELData[zonedata[z]->node_num];
368 ELAuxData* aux_new = new ELAuxData[pzonedata->node_num];
369 for(int i=0;i<pzonedata->node_num;i++)
370 {
371 fs_new[old_to_new[i]] = pzonedata->fs[i];
372 aux_new[old_to_new[i]] = pzonedata->aux[i];
373 }
374 delete [] pzonedata->fs;
375 delete [] pzonedata->aux;
376 pzonedata->fs = fs_new;
377 pzonedata->aux = aux_new;
378 }
379 if(zonedata[z]->material_type == Vacuum)
380 {
381 VacuumZone *pzonedata = dynamic_cast<VacuumZone *> (zonedata[z]);
382 VacuumData* fs_new = new VacuumData[zonedata[z]->node_num];
383 VacuumAuxData* aux_new = new VacuumAuxData[pzonedata->node_num];
384 for(int i=0;i<pzonedata->node_num;i++)
385 {
386 fs_new[old_to_new[i]] = pzonedata->fs[i];
387 aux_new[old_to_new[i]] = pzonedata->aux[i];
388 }
389 delete [] pzonedata->fs;
390 delete [] pzonedata->aux;
391 pzonedata->fs = fs_new;
392 pzonedata->aux = aux_new;
393 }
394 if(zonedata[z]->material_type == PML)
395 {
396 PMLZone *pzonedata = dynamic_cast<PMLZone *> (zonedata[z]);
397 PMLData* fs_new = new PMLData[zonedata[z]->node_num];
398 PMLAuxData* aux_new = new PMLAuxData[pzonedata->node_num];
399 for(int i=0;i<pzonedata->node_num;i++)
400 {
401 fs_new[old_to_new[i]] = pzonedata->fs[i];
402 aux_new[old_to_new[i]] = pzonedata->aux[i];
403 }
404 delete [] pzonedata->fs;
405 delete [] pzonedata->aux;
406 pzonedata->fs = fs_new;
407 pzonedata->aux = aux_new;
408 }
409 delete [] old_to_new;
410 }
411 build_zone();
412 }
413
414
415 int BSolver::build_least_squares()
416 {
417 double a=0,b=0,c=0,w=0;
418 for(int z=0; z<zone_num; z++)
419 if(zonedata[z]->material_type == Semiconductor)
420 {
421 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
422 for(int i=0; i<zone[z].davcell.size(); i++)
423 {
424 a=0,b=0,c=0,w=0;
425 VoronoiCell *pcell = zone[z].davcell.GetPointer(i);
426 for(int k=0;k<pcell->nb_num;k++)
427 {
428 const VoronoiCell *ncell = zone[z].davcell.GetPointer(pcell->nb_array[k]);
429 double dx = ncell->x - pcell->x;
430 double dy = ncell->y - pcell->y;
431 w=1.0/sqrt(dx*dx+dy*dy);
432 a+=w*w*dx*dx;
433 b+=w*w*dx*dy;
434 c+=w*w*dy*dy;
435 }
436 if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
437 {
438 InsulatorInterfaceBC *pbc;
439 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
440 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
441 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
442 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
443 const VoronoiCell* dcell = pz->pzone->davcell.GetPointer(n_node);
444 for(int k=1;k<dcell->nb_num-1;k++)
445 {
446 const VoronoiCell *ncell = pz->pzone->davcell.GetPointer(dcell->nb_array[k]);
447 double dx = ncell->x - dcell->x;
448 double dy = ncell->y - dcell->y;
449 w=1.0/sqrt(dx*dx+dy*dy);
450 a+=w*w*dx*dx;
451 b+=w*w*dx*dy;
452 c+=w*w*dy*dy;
453 }
454 }
455 pcell->sa = a/(a*c-b*b);
456 pcell->sb = b/(a*c-b*b);
457 pcell->sc = c/(a*c-b*b);
458 }
459 }
460 return 0;
461 }
462
463 int VacuumZone::Init(ZONE* zone, double Tl,PhysicalUnitScale *scale)
464 {
465 for(int i=0;i<node_num; i++)
466 {
467 aux[i].eps = mt->eps0;
468 aux[i].mu = mt->mu0;
469 }
470 return 0;
471 }
472
473 int PMLZone::Init(ZONE* zone, double Tl,PhysicalUnitScale *scale)
474 {
475 for(int i=0;i<node_num; i++)
476 {
477 aux[i].eps = mt->eps0;
478 aux[i].mu = mt->mu0;
479 }
480 return 0;
481 }
482
483 int ISZone::Init(ZONE* zone, double Tl,PhysicalUnitScale *scale)
484 {
485 for(int i=0;i<node_num; i++)
486 {
487 mt->mapping(&pzone->danode[i],&aux[i],0);
488 aux[i].Ex = aux[i].Ey = 0;
489 aux[i].affinity = mt->basic->Affinity(Tl);
490 aux[i].density = mt->basic->Density(Tl);
491 aux[i].eps = mt->eps0*mt->basic->Permittivity();
492 aux[i].mu = mt->mu0*mt->basic->Permeability();
493 fs[i].T = Tl;
494 fs[i].P = -aux[i].affinity;
495 }
496 return 0;
497 }
498
499 int ElZone::Init(ZONE* zone, double Tl,PhysicalUnitScale *scale)
500 {
501 for(int i=0;i<node_num; i++)
502 {
503 aux[i].affinity = mt->basic->Affinity(Tl);
504 aux[i].density = mt->basic->Density(Tl);
505 aux[i].eps = mt->eps0*mt->basic->Permittivity();
506 aux[i].mu = mt->mu0*mt->basic->Permeability();
507 fs[i].T = Tl;
508 fs[i].P = -aux[i].affinity;
509 }
510 return 0;
511 }
512
513
514 int SMCZone::Init(ZONE* zone, double Tl,PhysicalUnitScale *scale)
515 {
516 PetscScalar e = mt->e;
517 PetscScalar kb = mt->kb;
518 for(int i=0;i<node_num; i++)
519 {
520 fs[i].T = Tl;
521 fs[i].Tn=fs[i].Tp=fs[i].T;
522 mt->mapping(&pzone->danode[i],&aux[i],0);
523 double nie = mt->band->nie(fs[i].T);
524 double Nc = mt->band->Nc(fs[i].T);
525 if((aux[i].Nd-aux[i].Na)>0) //n-type
526 {
527 fs[i].n=((aux[i].Nd-aux[i].Na)+sqrt(pow(aux[i].Nd-aux[i].Na,2)+4*nie*nie))/2;
528 fs[i].p=nie*nie/fs[i].n;
529 }
530 else //p-type
531 {
532 fs[i].p=((aux[i].Na-aux[i].Nd)+sqrt(pow(aux[i].Na-aux[i].Nd,2)+4*nie*nie))/2;
533 fs[i].n=nie*nie/fs[i].p;
534 }
535 aux[i].affinity = mt->basic->Affinity(fs[i].T);
536 aux[i].density = mt->basic->Density(fs[i].T);
537 aux[i].eps = mt->eps0*mt->basic->Permittivity();
538 aux[i].mu = mt->mu0*mt->basic->Permeability();
539 aux[i].Eg = mt->band->Eg(fs[i].T)-mt->band->EgNarrow(fs[i].T);
540 aux[i].Nc = mt->band->Nc(fs[i].T);
541 aux[i].Nv = mt->band->Nv(fs[i].T);
542 fs[i].P = kb*fs[i].T/e*asinh((aux[i].Nd-aux[i].Na)/(2*nie)) - kb*fs[i].T/e*log(Nc/nie)
543 - aux[i].affinity - mt->band->EgNarrowToEc(fs[i].T);
544 fs[i].Eqc = -(e*fs[i].P + aux[i].affinity + mt->band->EgNarrowToEc(fs[i].T) + kb*fs[i].T*log(aux[i].Nc));
545 fs[i].Eqv = -(e*fs[i].P + aux[i].affinity - mt->band->EgNarrowToEv(fs[i].T) - kb*fs[i].T*log(aux[i].Nv)+ aux[i].Eg);
546 aux[i].mun = aux[i].mup = 0;
547 }
548 return 0;
549 }
550
551 void SMCZone::report()
552 {
553 gss_log.string_buf()<<"\nThis is a semiconductor Zone.\n";
554 gss_log.record();
555 }
556
557 void ISZone::report()
558 {
559 gss_log.string_buf()<<"\nThis is a insulator Zone.\n";
560 gss_log.record();
561 }
562
563 void ElZone::report()
564 {
565 gss_log.string_buf()<<"\nThis is a electrode.\n";
566 gss_log.record();
567 }
568
569 void VacuumZone::report()
570 {
571 gss_log.string_buf()<<"\nThis is vacuum area.\n";
572 gss_log.record();
573 }
574
575 void PMLZone::report()
576 {
577 gss_log.string_buf()<<"\nThis is PML.\n";
578 gss_log.record();
579 }
580
581
582
General-purpose Semiconductor Simulator