?  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 / ddm2e / ddm_nt2e.cc
blobcc4b0e8de01d68cda2861aa864f7f4356d54d292
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: Feb 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 "ddm_nt2e.h"
22 #include "private/kspimpl.h"
23 #include "private/snesimpl.h"
24 #include "src/snes/impls/tr/tr.h"
25 #include "log.h"
26
27
28 /* ----------------------------------------------------------------------------
29 * form_function_pn: This function setup DDM equation F(x)=0
30 */
31 void DDM_Solver_L2E::form_function_pn_2E(PetscScalar *x,PetscScalar *f)
32 {
33 // compute flux along triangle edges. semiconductor zone only
34 for(int z=0;z<zone_num;z++)
35 if(zonedata[z]->material_type==Semiconductor)
36 {
37 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
38 for(int i=0;i<zone[z].datri.size();i++)
39 {
40 Tri *ptri = &zone[z].datri[i];
41 pzonedata->F2E_Tri_ddm(ptri,x,f,zofs);
42 }
43 }
44
45 for(int z=0;z<zone_num;z++)
46 {
47 //semiconductor zone
48 if(zonedata[z]->material_type==Semiconductor)
49 {
50 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
51
52 // process electrode.
53 for(int i=0;i<pzonedata->electrode.size();i++)
54 {
55 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
56 pzonedata->F2E_om_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
57 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
58 pzonedata->F2E_stk_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
59 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
60 pzonedata->F2E_ins_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
61 }
62
63 for(int i=0;i<zone[z].davcell.size();i++)
64 {
65 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
66 if(!pcell->bc_index)
67 pzonedata->F2E_ddm_inner(i,x,f,ODE_formula,zofs);
68 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==NeumannBoundary)
69 pzonedata->F2E_ddm_neumannbc(i,x,f,ODE_formula,zofs,bc);
70 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
71 {
72 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
73 if(pbc->pinterface)//the ohmic bc is an electrode region
74 {
75 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
76 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
77 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
78 pzonedata->F2E_ddm_ombc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
79 }
80 else //the ohmic bc is a segment
81 pzonedata->F2E_ddm_ombc_segment(i,x,f,ODE_formula,zofs,bc);
82 }
83 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
84 {
85 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
86 if(pbc->pinterface)//the Schottky bc is an electrode region
87 {
88 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
89 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
90 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
91 pzonedata->F2E_ddm_stkbc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
92 }
93 else //the Schottky bc is a segment
94 pzonedata->F2E_ddm_stkbc_segment(i,x,f,ODE_formula,zofs,bc);
95 }
96 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
97 pzonedata->F2E_ddm_insulator_gate(i,x,f,ODE_formula,zofs,bc);
98 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
99 {
100 InsulatorInterfaceBC *pbc;
101 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
102 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
103 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
104 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
105 pzonedata->F2E_ddm_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
106 //pzonedata->F2_ddm_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
107 }
108 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HomoInterface)
109 {
110 HomoInterfaceBC *pbc;
111 pbc = dynamic_cast<HomoInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
112 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
113 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
114 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
115 pzonedata->F2E_ddm_homojunction(i,x,f,ODE_formula,zofs,bc,pz,n_node);
116 }
117 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HeteroInterface)
118 {
119 HeteroInterfaceBC *pbc;
120 pbc = dynamic_cast<HeteroInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
121 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
122 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
123 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
124 pzonedata->F2E_ddm_heterojunction(i,x,f,ODE_formula,zofs,bc,pz,n_node);
125 }
126 else //prevent some mistakes caused by inexact bc
127 pzonedata->F2E_ddm_inner(i,x,f,ODE_formula,zofs);
128 }
129 }
130
131 // Insulator zone
132 if(zonedata[z]->material_type==Insulator)
133 {
134 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
135 pzonedata->F2_efield_update(x,zofs,bc,zonedata);
136 for(int i=0;i<zone[z].davcell.size();i++)
137 {
138 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
139 if(!pcell->bc_index)
140 pzonedata->F2_ddm_inner(i,x,f,ODE_formula,zofs);
141 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==NeumannBoundary)
142 pzonedata->F2_ddm_neumannbc(i,x,f,ODE_formula,zofs,bc);
143 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
144 {
145 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
146 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
147 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
148 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
149 pzonedata->F2_ddm_chargebc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
150 }
151 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
152 {
153 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
154 if(pbc->pinterface)//the gate bc is an electrode region
155 {
156 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
157 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
158 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
159 pzonedata->F2_ddm_gatebc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
160 }
161 else //the gate bc is a segment
162 pzonedata->F2_ddm_gatebc_segment(i,x,f,ODE_formula,zofs,bc);
163 }
164 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
165 {
166 ElectrodeInsulatorBC *pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
167 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
168 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
169 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
170 pzonedata->F2_ddm_electrode_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
171 }
172 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
173 {
174 InsulatorInterfaceBC *pbc;
175 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
176 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
177 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
178 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
179 pzonedata->F2_ddm_semiconductor_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
180 }
181 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
182 {
183 InsulatorInsulatorBC *pbc;
184 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
185 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
186 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
187 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
188 pzonedata->F2_ddm_insulator_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
189 }
190 else //prevent some mistakes caused by inexact bc
191 pzonedata->F2_ddm_inner(i,x,f,ODE_formula,zofs);
192 }
193 for(int i=0;i<pzonedata->electrode.size();i++)
194 {
195 if(bc[pzonedata->electrode[i]].BCType==GateContact)
196 pzonedata->F2_gate_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
197 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
198 pzonedata->F2_charge_electrode(i,x,f,ODE_formula,zofs,bc);
199 }
200 }
201
202
203 // Electrode zone
204 if(zonedata[z]->material_type==Conductor)
205 {
206 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
207 for(int i=0;i<zone[z].davcell.size();i++)
208 {
209 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
210 if(!pcell->bc_index)
211 pzonedata->F2_ddm_inner(i,x,f,ODE_formula,zofs);
212 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==NeumannBoundary)
213 pzonedata->F2_ddm_neumannbc(i,x,f,ODE_formula,zofs,bc);
214 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
215 {
216 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
217 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
218 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
219 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
220 pzonedata->F2_ddm_chargebc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
221 }
222 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
223 {
224 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
225 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
226 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
227 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
228 pzonedata->F2_ddm_gatebc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
229 }
230 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
231 {
232 ElectrodeInsulatorBC *pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
233 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
234 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
235 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
236 pzonedata->F2_ddm_elec_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
237 }
238 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
239 {
240 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
241 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
242 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
243 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
244 pzonedata->F2_ddm_ombc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
245 }
246 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
247 {
248 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
249 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
250 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
251 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
252 pzonedata->F2_ddm_stkbc_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
253 }
254 else //prevent some mistakes caused by inexact bc
255 pzonedata->F2_ddm_inner(i,x,f,ODE_formula,zofs);
256 }
257 }
258 }
259 }
260
261
262 /* ----------------------------------------------------------------------------
263 * form_jacobian_pn: This function setup jacobian matrix F'(x) of DDM equation F(x)=0
264 * dual carrier edition
265 */
266 void DDM_Solver_L2E::form_jacobian_pn_2E(PetscScalar *x, Mat *jac, Mat *jtmp)
267 {
268
269 //Compute Jacobian entries for flux along triangle edges.
270 for(int z=0;z<zone_num;z++)
271 {
272 if(zonedata[z]->material_type==Semiconductor)
273 {
274 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
275
276 // compute flux Jacobian along triangle edges
277 for(int i=0;i<zone[z].datri.size();i++)
278 {
279 Tri *ptri = &zone[z].datri[i];
280 pzonedata->J2E_Tri_ddm(ptri,x,jtmp,zofs);
281 }
282 }
283 }
284 MatAssemblyBegin(*jtmp,MAT_FINAL_ASSEMBLY);
285 MatAssemblyEnd(*jtmp,MAT_FINAL_ASSEMBLY);
286
287
288 //Compute Jacobian entries and insert into matrix.
289 for(int z=0;z<zone_num;z++)
290 {
291 if(zonedata[z]->material_type==Semiconductor)
292 {
293 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
294
295 // process electrode.
296 for(int i=0;i<pzonedata->electrode.size();i++)
297 {
298 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
299 pzonedata->J2E_om_electrode(i,x,jac,jtmp,ODE_formula,zofs,bc,DeviceDepth);
300 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
301 pzonedata->J2E_stk_electrode(i,x,jac,jtmp,ODE_formula,zofs,bc,DeviceDepth);
302 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
303 pzonedata->J2E_ins_electrode(i,x,jac,jtmp,ODE_formula,zofs,bc,DeviceDepth);
304 }
305
306 for(int i=0;i<zone[z].davcell.size();i++)
307 {
308 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
309 SemiData *pfs = &pzonedata->fs[i];
310 if(!pcell->bc_index)
311 pzonedata->J2E_ddm_inner(i,x,jac,jtmp,ODE_formula,zofs);
312 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==NeumannBoundary)
313 pzonedata->J2E_ddm_neumannbc(i,x,jac,jtmp,ODE_formula,zofs,bc);
314 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
315 {
316 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
317 if(pbc->pinterface)//the ohmic bc is an electrode region
318 {
319 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
320 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
321 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
322 pzonedata->J2E_ddm_ombc_interface(i,x,jac,jtmp,ODE_formula,zofs,bc,pz,n_node);
323 }
324 else //the ohmic bc is a segment
325 pzonedata->J2E_ddm_ombc_segment(i,x,jac,jtmp,ODE_formula,zofs,bc);
326 }
327 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
328 {
329 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
330 if(pbc->pinterface)//the Schottky bc is an electrode region
331 {
332 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
333 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
334 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
335 pzonedata->J2E_ddm_stkbc_interface(i,x,jac,jtmp,ODE_formula,zofs,bc,pz,n_node);
336 }
337 else //the Schottky bc is a segment
338 pzonedata->J2E_ddm_stkbc_segment(i,x,jac,jtmp,ODE_formula,zofs,bc);
339 }
340 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
341 pzonedata->J2E_ddm_insulator_gate(i,x,jac,jtmp,ODE_formula,zofs,bc);
342 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
343 {
344 InsulatorInterfaceBC *pbc;
345 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
346 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
347 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
348 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
349 pzonedata->J2E_ddm_interface(i,x,jac,jtmp,ODE_formula,zofs,bc,pz,n_node);
350 //pzonedata->J2_ddm_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
351 }
352 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HomoInterface)
353 {
354 HomoInterfaceBC *pbc;
355 pbc = dynamic_cast<HomoInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
356 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
357 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
358 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
359 pzonedata->J2E_ddm_homojunction(i,x,jac,jtmp,ODE_formula,zofs,bc,pz,n_node);
360 }
361 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HeteroInterface)
362 {
363 HeteroInterfaceBC *pbc;
364 pbc = dynamic_cast<HeteroInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
365 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
366 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
367 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
368 pzonedata->J2E_ddm_heterojunction(i,x,jac,jtmp,ODE_formula,zofs,bc,pz,n_node);
369 }
370 else //prevent some mistakes caused by inexact bc
371 pzonedata->J2E_ddm_inner(i,x,jac,jtmp,ODE_formula,zofs);
372 }
373 }
374
375 // Insulator zone
376 if(zonedata[z]->material_type==Insulator)
377 {
378 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
379 for(int i=0;i<zone[z].davcell.size();i++)
380 {
381 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
382 if(!pcell->bc_index)
383 pzonedata->J2_ddm_inner(i,x,jac,ODE_formula,zofs);
384 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==NeumannBoundary)
385 pzonedata->J2_ddm_neumannbc(i,x,jac,ODE_formula,zofs,bc);
386 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
387 {
388 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
389 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
390 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
391 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
392 pzonedata->J2_ddm_chargebc_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
393 }
394 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
395 {
396 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
397 if(pbc->pinterface)//the gate bc is an electrode region
398 {
399 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
400 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
401 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
402 pzonedata->J2_ddm_gatebc_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
403 }
404 else //the gate bc is a segment
405 pzonedata->J2_ddm_gatebc_segment(i,x,jac,ODE_formula,zofs,bc);
406 }
407 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
408 {
409 ElectrodeInsulatorBC *pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
410 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
411 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
412 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
413 pzonedata->J2_ddm_electrode_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
414 }
415 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
416 {
417 InsulatorInterfaceBC *pbc;
418 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
419 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
420 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
421 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
422 pzonedata->J2_ddm_semiconductor_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
423 }
424 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
425 {
426 InsulatorInsulatorBC *pbc;
427 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
428 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
429 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
430 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
431 pzonedata->J2_ddm_insulator_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
432 }
433 else //prevent some mistakes caused by inexact bc
434 pzonedata->J2_ddm_inner(i,x,jac,ODE_formula,zofs);
435 }
436 for(int i=0;i<pzonedata->electrode.size();i++)
437 {
438 if(bc[pzonedata->electrode[i]].BCType==GateContact)
439 pzonedata->J2_gate_electrode(i,x,jac,ODE_formula,zofs,bc,DeviceDepth);
440 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
441 pzonedata->J2_charge_electrode(i,x,jac,ODE_formula,zofs,bc);
442 }
443 }
444
445 // Electrode zone
446 if(zonedata[z]->material_type==Conductor)
447 {
448 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
449 for(int i=0;i<zone[z].davcell.size();i++)
450 {
451 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
452 if(!pcell->bc_index)
453 pzonedata->J2_ddm_inner(i,x,jac,ODE_formula,zofs);
454 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==NeumannBoundary)
455 pzonedata->J2_ddm_neumannbc(i,x,jac,ODE_formula,zofs,bc);
456 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
457 {
458 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
459 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
460 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
461 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
462 pzonedata->J2_ddm_chargebc_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
463 }
464 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
465 {
466 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
467 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
468 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
469 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
470 pzonedata->J2_ddm_gatebc_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
471 }
472 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
473 {
474 ElectrodeInsulatorBC *pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
475 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
476 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
477 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
478 pzonedata->J2_ddm_elec_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
479 }
480 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
481 {
482 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
483 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
484 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
485 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
486 pzonedata->J2_ddm_ombc_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
487 }
488 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
489 {
490 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
491 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
492 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
493 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
494 pzonedata->J2_ddm_stkbc_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
495 }
496 else //prevent some mistakes caused by inexact bc
497 pzonedata->J2_ddm_inner(i,x,jac,ODE_formula,zofs);
498 }
499 }
500
501 }
502 }
503
504
505
506 /* ----------------------------------------------------------------------------
507 * SNES_form_function_pn_2E: wrap function for petsc nonlinear solver
508 */
509 PetscErrorCode SNES_form_function_pn_2E(SNES snes, Vec x,Vec f,void *dummy)
510 {
511 PetscScalar *xx,*ff;
512 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)dummy;
513 VecZeroEntries(f);
514 //Get pointers to vector data.
515 VecGetArray(x,&xx);
516 VecGetArray(f,&ff);
517
518 ps->form_function_pn_2E(xx,ff);
519 ps->error_norm_estimat(xx,ff);
520
521 //Restore vectors
522 VecRestoreArray(x,&xx);
523 VecRestoreArray(f,&ff);
524 VecNorm(f,NORM_2,&ps->norm);
525
526 return 0;
527 }
528
529
530 /* ----------------------------------------------------------------------------
531 * SNES_form_jacobian_pn_2E: wrap function for petsc nonlinear solver
532 */
533 PetscErrorCode SNES_form_jacobian_pn_2E(SNES snes, Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
534 {
535 PetscScalar *xx;
536 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)dummy;
537 //Get pointer to vector data
538 VecGetArray(x,&xx);
539 //clear old matrix
540 MatZeroEntries(*jac);
541 MatZeroEntries(ps->JTmp);
542 //build matrix here
543 ps->form_jacobian_pn_2E(xx,jac,&ps->JTmp);
544 *flag = SAME_NONZERO_PATTERN;
545 //Restore vector
546 VecRestoreArray(x,&xx);
547 //Assemble matrix
548 MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
549 MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
550 //MatView(*jac,PETSC_VIEWER_DRAW_WORLD);
551 //getchar();
552 return 0;
553 }
554
555
556 /* ----------------------------------------------------------------------------
557 * SNESMF_form_jacobian_pn_2E: wrap function for petsc nonlinear solver with
558 * matrix-free method
559 */
560 PetscErrorCode SNESMF_form_jacobian_pn_2E(SNES snes, Vec x,Mat *A,Mat *B,MatStructure *flag,void *dummy)
561 {
562 PetscScalar *xx;
563 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)dummy;
564 //Get pointer to vector data
565 VecGetArray(x,&xx);
566 //clear old matrix
567 MatZeroEntries(*B);
568 MatZeroEntries(ps->JTmp);
569 //build matrix here
570 ps->form_jacobian_pn_2E(xx,B,&ps->JTmp);
571 *flag = SAME_NONZERO_PATTERN;
572 //Restore vector
573 VecRestoreArray(x,&xx);
574 //Assemble matrix
575 MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
576 MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
577 MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
578 MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
579 //MatView(*jac,PETSC_VIEWER_DRAW_WORLD);
580 //getchar();
581 return 0;
582 }
583
584
585 /* ----------------------------------------------------------------------------
586 * error_norm_pn_2E: This function compute X and RHS error norms.
587 */
588 void DDM_Solver_L2E:: error_norm_estimat(PetscScalar *x,PetscScalar *f)
589 {
590 // do clear
591 potential_norm = 0;
592 electron_norm = 0;
593 hole_norm = 0;
594 temperature_norm = 0;
595 possion_norm = 0;
596 elec_continuty_norm = 0;
597 hole_continuty_norm = 0;
598 electrode_norm = 0;
599 heat_equation_norm = 0;
600
601 int offset=0;
602 for(int z=0;z<zone_num;z++)
603 {
604 if(zonedata[z]->material_type==Semiconductor)
605 {
606 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
607 for(int i=0;i<zone[z].davcell.size();i++)
608 {
609 potential_norm += x[offset+0]*x[offset+0];
610 electron_norm += x[offset+1]*x[offset+1];
611 hole_norm += x[offset+2]*x[offset+2];
612 temperature_norm += x[offset+3]*x[offset+3];
613 possion_norm += f[offset+0]*f[offset+0];
614 elec_continuty_norm += f[offset+1]*f[offset+1];
615 hole_continuty_norm += f[offset+2]*f[offset+2];
616 heat_equation_norm += f[offset+3]*f[offset+3];
617 offset += 4;
618 }
619 for(int i=0;i<pzonedata->electrode.size();i++)
620 {
621 potential_norm += x[offset]*x[offset];
622 electrode_norm += f[offset]*f[offset];
623 offset += 1;
624 }
625 }
626 if(zonedata[z]->material_type==Insulator)
627 {
628 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
629 for(int i=0;i<zone[z].davcell.size();i++)
630 {
631 potential_norm += x[offset]*x[offset];
632 temperature_norm += x[offset+1]*x[offset+1];
633 possion_norm += f[offset]*f[offset];
634 heat_equation_norm += f[offset+1]*f[offset+1];
635 offset += 2;
636 }
637 for(int i=0;i<pzonedata->electrode.size();i++)
638 {
639 potential_norm += x[offset]*x[offset];
640 electrode_norm += f[offset]*f[offset];
641 offset += 1;
642 }
643 }
644 if(zonedata[z]->material_type==Conductor)
645 {
646 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
647 for(int i=0;i<zone[z].davcell.size();i++)
648 {
649 potential_norm += x[offset]*x[offset];
650 temperature_norm += x[offset+1]*x[offset+1];
651 possion_norm += f[offset]*f[offset];
652 heat_equation_norm += f[offset+1]*f[offset+1];
653 offset += 2;
654 }
655 }
656 }
657
658 potential_norm = sqrt(potential_norm);
659 electron_norm = sqrt(electron_norm);
660 hole_norm = sqrt(hole_norm);
661 temperature_norm = sqrt(temperature_norm);
662 possion_norm = sqrt(possion_norm);
663 elec_continuty_norm = sqrt(elec_continuty_norm);
664 hole_continuty_norm = sqrt(hole_continuty_norm);
665 electrode_norm = sqrt(electrode_norm);
666 heat_equation_norm = sqrt(heat_equation_norm);
667
668 }
669
670
671 inline PetscErrorCode LineSearch_ConvergenceTest_L2E(SNES snes,PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm,
672 SNESConvergedReason *reason,void *dummy)
673 {
674 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)dummy;
675
676 *reason = SNES_CONVERGED_ITERATING;
677 if (!it)
678 {
679 snes->ttol = fnorm*snes->rtol;
680 gss_log.string_buf()<<" "<<"its\t"<<"| Eq(V) | "<<"| Eq(n) | "<<"| Eq(p) | "<<"| Eq(T) | "<<"|delta x|\n"
681 <<"----------------------------------------------------------------------\n";
682 gss_log.record();
683 }
684
685 gss_log.string_buf().precision(3);
686 gss_log.string_buf()<<" "<<it<<"\t"
687 <<ps->possion_norm<<" "
688 <<ps->elec_continuty_norm<<" "
689 <<ps->hole_continuty_norm<<" "
690 <<ps->heat_equation_norm<<" "
691 <<pnorm<<"\n";
692 gss_log.record();
693 gss_log.string_buf().precision(6);
694
695 if (fnorm != fnorm)
696 {
697 *reason = SNES_DIVERGED_FNORM_NAN;
698 }
699 else if (ps->possion_norm < ps->possion_abs_toler &&
700 ps->elec_continuty_norm < ps->elec_continuty_abs_toler &&
701 ps->hole_continuty_norm < ps->hole_continuty_abs_toler &&
702 ps->electrode_norm < ps->electrode_abs_toler &&
703 ps->heat_equation_norm < ps->heat_equation_abs_toler)
704 {
705 *reason = SNES_CONVERGED_FNORM_ABS;
706 }
707 else if (snes->nfuncs >= snes->max_funcs)
708 {
709 *reason = SNES_DIVERGED_FUNCTION_COUNT;
710 }
711
712 if (it && !*reason)
713 {
714 if (fnorm <= snes->ttol &&
715 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
716 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
717 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
718 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler &&
719 ps->heat_equation_norm < ps->toler_relax*ps->heat_equation_abs_toler)
720 {
721 *reason = SNES_CONVERGED_FNORM_RELATIVE;
722 }
723 else if (pnorm < ps->relative_toler &&
724 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
725 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
726 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
727 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler &&
728 ps->heat_equation_norm < ps->toler_relax*ps->heat_equation_abs_toler)
729 {
730 *reason = SNES_CONVERGED_PNORM_RELATIVE;
731 }
732 }
733
734 if(*reason>0)
735 {
736 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
737 <<" "<<SNESConvergedReasons[*reason]<<" *residual norm = "<<fnorm<<"\n\n\n";
738 gss_log.record();
739 }
740
741 return(0);
742 }
743
744
745 inline PetscErrorCode TrustRegion_ConvergenceTest_L2E(SNES snes,PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm,
746 SNESConvergedReason *reason,void *dummy)
747 {
748 SNES_TR *neP = (SNES_TR *)snes->data;
749 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)dummy;
750
751 if (!ps->its)
752 {
753 gss_log.string_buf()<<" "<<"its\t"<<"| Eq(V) | "<<"| Eq(n) | "<<"| Eq(p) | "<<"| Eq(T) | "<<"|delta x|\n"
754 <<"----------------------------------------------------------------------\n";
755 gss_log.record();
756 }
757
758 gss_log.string_buf().precision(3);
759 gss_log.string_buf()<<" "<<ps->its++<<"\t"
760 <<ps->possion_norm<<" "
761 <<ps->elec_continuty_norm<<" "
762 <<ps->hole_continuty_norm<<" "
763 <<ps->heat_equation_norm<<" "
764 <<pnorm<<"\n";
765 gss_log.record();
766 gss_log.string_buf().precision(6);
767
768 *reason = SNES_CONVERGED_ITERATING;
769
770 if (!it)
771 {
772 snes->ttol = fnorm*snes->rtol;
773 }
774
775 if (fnorm != fnorm)
776 {
777 *reason = SNES_DIVERGED_FNORM_NAN;
778 }
779 else if (neP->delta < xnorm * snes->deltatol )
780 {
781 if( ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
782 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
783 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
784 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler &&
785 ps->heat_equation_norm < ps->toler_relax*ps->heat_equation_abs_toler)
786 {
787 *reason = SNES_CONVERGED_TR_DELTA;
788 }
789 else
790 {
791 *reason = SNES_DIVERGED_LOCAL_MIN;
792 }
793 }
794 else if (ps->possion_norm < ps->possion_abs_toler &&
795 ps->elec_continuty_norm < ps->elec_continuty_abs_toler &&
796 ps->hole_continuty_norm < ps->hole_continuty_abs_toler &&
797 ps->electrode_norm < ps->electrode_abs_toler &&
798 ps->heat_equation_norm < ps->heat_equation_abs_toler)
799 {
800 *reason = SNES_CONVERGED_FNORM_ABS;
801 }
802 else if (snes->nfuncs >= snes->max_funcs)
803 {
804 *reason = SNES_DIVERGED_FUNCTION_COUNT;
805 }
806
807 if (it && !*reason)
808 {
809 if (fnorm <= snes->ttol)
810 {
811 *reason = SNES_CONVERGED_FNORM_RELATIVE;
812 }
813 else if (pnorm < ps->relative_toler &&
814 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
815 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
816 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
817 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler &&
818 ps->heat_equation_norm < ps->toler_relax*ps->heat_equation_abs_toler)
819 {
820 *reason = SNES_CONVERGED_PNORM_RELATIVE;
821 }
822 }
823
824 if (snes->nfuncs >= snes->max_funcs)
825 {
826 *reason = SNES_DIVERGED_FUNCTION_COUNT;
827 }
828 if(*reason>0)
829 {
830 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
831 <<" "<<SNESConvergedReasons[*reason]<<" *residual norm = "<<fnorm<<"\n\n\n";
832 gss_log.record();
833 }
834
835 return(0);
836 }
837
838
839 PetscErrorCode LimitorByBankRose_L2E(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
840 {
841 int it;
842 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)checkctx;
843 PetscScalar *xx;
844 PetscScalar *yy;
845 PetscScalar *ww;
846 VecGetArray(x,&xx);
847 VecGetArray(y,&yy);
848 VecGetArray(w,&ww);
849
850 SNESGetIterationNumber(snes,&it);
851 static PetscScalar K;
852 if(it==0) K=0;
853
854 Vec gk0,gk1;
855 int flag=0;
856 VecDuplicate(x,&gk0);
857 VecDuplicate(x,&gk1);
858 SNESComputeFunction(snes,x,gk0);
859 PetscScalar gk0_norm;
860 VecNorm(gk0,NORM_2,&gk0_norm);
861 while(1)
862 {
863 PetscScalar tk=1.0/(1+K*gk0_norm);
864 if(tk<1e-4) break;
865 //w=x-tk*y;
866 flag=1;
867 VecWAXPY(w,-tk,y,x);
868 SNESComputeFunction(snes,w,gk1);
869 PetscScalar gk1_norm;
870 VecNorm(gk1,NORM_2,&gk1_norm);
871 if((1-gk1_norm/gk0_norm)<0.9*tk)
872 {
873 if(K==0) K=1;
874 else K*=10;
875 }
876 else
877 {
878 K/=10;
879 break;
880 }
881 }
882
883 //prevent negative carrier density and temperature underflow
884 int offset=0;
885 for(int z=0;z<ps->zone_num;z++)
886 {
887 if(ps->zonedata[z]->material_type==Semiconductor)
888 {
889 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
890 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
891 {
892 if(fabs(yy[offset])>1.0) {ww[offset] = xx[offset] - dsign(yy[offset])*1.0;flag=1;}
893 if (ww[offset+1] <0 ) {ww[offset+1]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
894 if (ww[offset+2] <0 ) {ww[offset+2]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
895 if (ww[offset+3] <0.7*ps->LatticeTemp ) {ww[offset+3]=0.7*ps->LatticeTemp;flag=1;}
896 offset += 4;
897 }
898 offset += pzonedata->electrode.size();
899 }
900 else if(ps->zonedata[z]->material_type==Insulator)
901 {
902 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
903 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
904 {
905 if (ww[offset+1]<0.7*ps->LatticeTemp ) {ww[offset+1]=0.7*ps->LatticeTemp;flag=1;}
906 offset += 2;
907 }
908 offset += pzonedata->electrode.size();
909 }
910 else if(ps->zonedata[z]->material_type==Conductor)
911 {
912 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
913 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
914 {
915 if (ww[offset+1]<0.7*ps->LatticeTemp ) {ww[offset+1]=0.7*ps->LatticeTemp;flag=1;}
916 offset += 2;
917 }
918 }
919 }
920
921 if(flag)
922 {
923 *changed_y = PETSC_FALSE;
924 *changed_w = PETSC_TRUE;
925 }
926 VecDestroy(gk0);
927 VecDestroy(gk1);
928 VecRestoreArray(x,&xx);
929 VecRestoreArray(y,&yy);
930 VecRestoreArray(w,&ww);
931 return(0);
932 }
933
934 PetscErrorCode LimitorByPotential_L2E(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
935 {
936 PetscScalar *xx;
937 PetscScalar *yy;
938 PetscScalar *ww;
939 VecGetArray(x,&xx);
940 VecGetArray(y,&yy);
941 VecGetArray(w,&ww);
942 PetscScalar dV_max=0;
943
944 //search for dV_max;
945 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)checkctx;
946 int offset=0;
947 for(int z=0;z<ps->zone_num;z++)
948 {
949 if(ps->zonedata[z]->material_type==Semiconductor)
950 {
951 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
952 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
953 {
954 if(fabs(yy[offset])>dV_max) {dV_max=fabs(yy[offset]); }
955 offset += 4;
956 }
957 offset += pzonedata->electrode.size();
958 }
959 else if(ps->zonedata[z]->material_type==Insulator)
960 {
961 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
962 offset += 2*pzonedata->pzone->davcell.size();
963 offset += pzonedata->electrode.size();
964 }
965 else if(ps->zonedata[z]->material_type==Conductor)
966 {
967 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
968 offset += 2*pzonedata->pzone->davcell.size();
969 }
970 }
971 if(dV_max<1e-2) return(0);
972
973 PetscScalar Vt=0.026*ps->LatticeTemp;
974 PetscScalar f = log(1+dV_max/Vt)/(dV_max/Vt);
975
976 //prevent negative carrier density and temperature underflow
977 offset=0;
978 for(int z=0;z<ps->zone_num;z++)
979 {
980 if(ps->zonedata[z]->material_type==Semiconductor)
981 {
982 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
983 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
984 {
985 ww[offset] = xx[offset] - f*yy[offset];
986 if (ww[offset+1] <0 ) ww[offset+1]=1.0*pow(ps->scale_unit.s_centimeter,-3);
987 if (ww[offset+2] <0 ) ww[offset+2]=1.0*pow(ps->scale_unit.s_centimeter,-3);
988 if (ww[offset+3] <0.7*ps->LatticeTemp ) ww[offset+3]=0.7*ps->LatticeTemp;
989 offset += 4;
990 }
991 offset += pzonedata->electrode.size();
992 }
993 else if(ps->zonedata[z]->material_type==Insulator)
994 {
995 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
996 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
997 {
998 ww[offset] = xx[offset] - f*yy[offset];
999 if (ww[offset+1]<0.7*ps->LatticeTemp ) ww[offset+1]=0.7*ps->LatticeTemp;
1000 offset += 2;
1001 }
1002 offset += pzonedata->electrode.size();
1003 }
1004 else if(ps->zonedata[z]->material_type==Conductor)
1005 {
1006 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
1007 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
1008 {
1009 ww[offset] = xx[offset] - f*yy[offset];
1010 if (ww[offset+1]<0.7*ps->LatticeTemp ) ww[offset+1]=0.7*ps->LatticeTemp;
1011 offset += 2;
1012 }
1013 }
1014 }
1015
1016 VecRestoreArray(x,&xx);
1017 VecRestoreArray(y,&yy);
1018 VecRestoreArray(w,&ww);
1019 *changed_y = PETSC_FALSE;
1020 *changed_w = PETSC_TRUE;
1021 return(0);
1022 }
1023
1024
1025 PetscErrorCode LimitorNonNegativeCarrier_L2E(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
1026 {
1027 PetscScalar *xx;
1028 PetscScalar *yy;
1029 PetscScalar *ww;
1030 VecGetArray(x,&xx);
1031 VecGetArray(y,&yy);
1032 VecGetArray(w,&ww);
1033 PetscScalar WorstRatio=0.0;
1034
1035 DDM_Solver_L2E *ps = (DDM_Solver_L2E*)checkctx;
1036 int flag=0;
1037 int offset=0;
1038 for(int z=0;z<ps->zone_num;z++)
1039 {
1040 if(ps->zonedata[z]->material_type==Semiconductor)
1041 {
1042 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
1043 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
1044 {
1045 PetscScalar r;
1046 if(fabs(yy[offset])>1.0) {ww[offset] = xx[offset] - dsign(yy[offset])*1.0;flag=1;}
1047 if (ww[offset+1]<0.0) {ww[offset+1]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
1048 if (ww[offset+2]<0.0) {ww[offset+2]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
1049 if (ww[offset+3]<0.7*ps->LatticeTemp ) {ww[offset+3]=0.7*ps->LatticeTemp;flag=1;}
1050 offset += 4;
1051 }
1052 offset += pzonedata->electrode.size();
1053 }
1054 else if(ps->zonedata[z]->material_type==Insulator)
1055 {
1056 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
1057 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
1058 {
1059 if (ww[offset+1]<0.7*ps->LatticeTemp ) {ww[offset+1]=0.7*ps->LatticeTemp;flag=1;}
1060 offset += 2;
1061 }
1062 offset += pzonedata->electrode.size();
1063 }
1064 else if(ps->zonedata[z]->material_type==Conductor)
1065 {
1066 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
1067 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
1068 {
1069 if (ww[offset+1]<0.7*ps->LatticeTemp ) {ww[offset+1]=0.7*ps->LatticeTemp;flag=1;}
1070 offset += 2;
1071 }
1072 }
1073 }
1074 VecRestoreArray(x,&xx);
1075 VecRestoreArray(y,&yy);
1076 VecRestoreArray(w,&ww);
1077
1078 if(flag)
1079 {
1080 *changed_y = PETSC_FALSE;
1081 *changed_w = PETSC_TRUE;
1082 }
1083 return(0);
1084 }
1085
1086
1087 /* ----------------------------------------------------------------------------
1088 * DDM_Solver_L2E::init_solver: This function do initial setup for nonlinear solver
1089 */
1090 int DDM_Solver_L2E::init_solver(SolveDefine &sv)
1091 {
1092 gss_log.string_buf()<<"DDM solver Level 2E init...";
1093 //set Tolerances
1094 relative_toler = sv.relative_toler;
1095 toler_relax = sv.toler_relax;
1096 possion_abs_toler = sv.possion_abs_toler;
1097 elec_continuty_abs_toler = sv.elec_continuty_abs_toler;
1098 hole_continuty_abs_toler = sv.hole_continuty_abs_toler;
1099 heat_equation_abs_toler = sv.heat_equation_abs_toler;
1100 electrode_abs_toler = sv.electrode_abs_toler;
1101 //if user defined a matrix-free method
1102 PetscOptionsHasName(PETSC_NULL,"-snes_mf_operator",&mf_flg);
1103
1104 // compute the scale of problem
1105 zofs.resize(zone_num+1);
1106 zofs[0] = 0;
1107 for(int i=0;i<zone_num;i++)
1108 {
1109 if(zonedata[i]->material_type==Semiconductor) //semiconductor zone
1110 {
1111 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1112 N+= 4*pzonedata->pzone->davcell.size();
1113 N+=pzonedata->electrode.size(); //additional equation for electrode
1114 zofs[i+1] = N;
1115 }
1116 else if(zonedata[i]->material_type==Insulator) //Insulator zone
1117 {
1118 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1119 N += 2*zone[i].davcell.size();
1120 N+=pzonedata->electrode.size(); //additional equation for electrode
1121 zofs[i+1] = N;
1122 }
1123 else if(zonedata[i]->material_type==Conductor) //Electrode zone
1124 {
1125 N += 2*zone[i].davcell.size();
1126 zofs[i+1] = N;
1127 }
1128 else //other zones, such as PML and vacuum
1129 {
1130 zofs[i+1] = N;
1131 }
1132 }
1133
1134 VecCreateSeq(PETSC_COMM_SELF,N,&x);
1135 VecDuplicate(x,&r);
1136 VecDuplicate(x,&x_n);
1137 VecDuplicate(x,&x_n1);
1138 VecDuplicate(x,&x_n2);
1139 VecDuplicate(x,&xp);
1140 VecDuplicate(x,&LTE);
1141
1142 //Evaluate initial guess,importnat for newton solver
1143 PetscScalar init_value[4];
1144 int index[4];
1145 int offset=0;
1146 for(int z=0;z<zone_num;z++)
1147 {
1148 if(zonedata[z]->material_type==Semiconductor)
1149 {
1150 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1151 for(int i=0;i<zone[z].davcell.size();i++)
1152 {
1153 index[0] = offset+0;
1154 index[1] = offset+1;
1155 index[2] = offset+2;
1156 index[3] = offset+3;
1157 init_value[0] = pzonedata->fs[i].P;
1158 init_value[1] = pzonedata->fs[i].n;
1159 init_value[2] = pzonedata->fs[i].p;
1160 init_value[3] = pzonedata->fs[i].T;
1161 VecSetValues(x,4,index,init_value,INSERT_VALUES);
1162 offset += 4;
1163 }
1164 for(int j=0;j<pzonedata->electrode.size();j++)
1165 {
1166 VecSetValue(x,offset,bc.Get_pointer(pzonedata->electrode[j])->Get_Potential(),INSERT_VALUES);
1167 offset += 1;
1168 }
1169 }
1170 if(zonedata[z]->material_type==Insulator)
1171 {
1172 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1173 for(int i=0;i<zone[z].davcell.size();i++)
1174 {
1175 VecSetValue(x,offset++,pzonedata->fs[i].P,INSERT_VALUES);
1176 VecSetValue(x,offset++,pzonedata->fs[i].T,INSERT_VALUES);
1177 }
1178 for(int j=0;j<pzonedata->electrode.size();j++)
1179 {
1180 VecSetValue(x,offset,bc.Get_pointer(pzonedata->electrode[j])->Get_Potential(),INSERT_VALUES);
1181 offset += 1;
1182 }
1183 }
1184 if(zonedata[z]->material_type==Conductor)
1185 {
1186 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1187 for(int i=0;i<zone[z].davcell.size();i++)
1188 {
1189 VecSetValue(x,offset++,pzonedata->fs[i].P,INSERT_VALUES);
1190 VecSetValue(x,offset++,pzonedata->fs[i].T,INSERT_VALUES);
1191 }
1192 }
1193 }
1194 //
1195 VecAssemblyBegin(x);
1196 VecAssemblyEnd(x);
1197
1198 // Create Jacobian matrix data structure
1199 // pre-alloc approximate memory
1200 int *nnz = new int[N];
1201 int *p = nnz;
1202 for(int i=0;i<zone_num;i++)
1203 {
1204 //semiconductor zone
1205 if(zonedata[i]->material_type==Semiconductor)
1206 {
1207 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1208 for(int j=0;j<pzonedata->pzone->davcell.size();j++)
1209 {
1210 const VoronoiCell* pcell = pzonedata->pzone->davcell.GetPointer(j);
1211 //inner node, alloc exact memory.
1212 if(!pcell->bc_index || bc[pcell->bc_index-1].psegment->interface==-1)
1213 {
1214 *p++ = 4*(pzonedata->pzone->davcell[j].nb_num+1);
1215 *p++ = 4*(pzonedata->pzone->davcell[j].nb_num+1);
1216 *p++ = 4*(pzonedata->pzone->davcell[j].nb_num+1);
1217 *p++ = 4*(pzonedata->pzone->davcell[j].nb_num+1);
1218 }
1219 else //interface node, slightly more than needed.
1220 {
1221 *p++ = 4*(2*pzonedata->pzone->davcell[j].nb_num+4);
1222 *p++ = 4*(2*pzonedata->pzone->davcell[j].nb_num+4);
1223 *p++ = 4*(2*pzonedata->pzone->davcell[j].nb_num+4);
1224 *p++ = 4*(2*pzonedata->pzone->davcell[j].nb_num+4);
1225 }
1226 }
1227 for(int j=0;j<pzonedata->electrode.size();j++)
1228 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
1229 }
1230 //Insulator zones, poisson equation and heat transfer equation
1231 if(zonedata[i]->material_type==Insulator)
1232 {
1233 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1234 for(int j=0;j<zone[i].davcell.size();j++)
1235 {
1236 *p++ = 2*(zone[i].davcell[j].nb_num+1);
1237 *p++ = 2*(zone[i].davcell[j].nb_num+1);
1238 }
1239 for(int j=0;j<pzonedata->electrode.size();j++)
1240 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
1241 }
1242 //Electrode zones, poisson equation and heat transfer equation
1243 if(zonedata[i]->material_type==Conductor)
1244 {
1245 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
1246 for(int j=0;j<zone[i].davcell.size();j++)
1247 {
1248 *p++ = 2*(zone[i].davcell[j].nb_num+1);
1249 *p++ = 2*(zone[i].davcell[j].nb_num+1);
1250 }
1251 }
1252 }
1253
1254 MatCreate(PETSC_COMM_SELF,&J);
1255 MatSetSizes(J,N,N,N,N);
1256 if(sv.LS=="superlu")
1257 {
1258 #ifdef PETSC_HAVE_SUPERLU
1259 MatSetType(J,MATSUPERLU);
1260 #else
1261 MatSetType(J,MATSEQAIJ);
1262 #endif
1263 }
1264 else if(sv.LS=="umfpack")
1265 {
1266 #ifdef PETSC_HAVE_UMFPACK
1267 MatSetType(J,MATUMFPACK);
1268 #else
1269 MatSetType(J,MATSEQAIJ);
1270 #endif
1271 }
1272 else
1273 {
1274 MatSetType(J,MATSEQAIJ);
1275 }
1276 MatSeqAIJSetPreallocation(J,0,nnz);
1277
1278 MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,0,nnz,&JTmp);
1279 if(mf_flg)
1280 MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,0,nnz,&JPrec);
1281 delete [] nnz;
1282
1283
1284 SNESCreate(PETSC_COMM_WORLD,&snes);
1285 SNESSetFunction(snes,r,SNES_form_function_pn_2E,this);
1286
1287 //if user defined a matrix-free method
1288 if(mf_flg)
1289 {
1290 gss_log.string_buf()<<"Matrix-Free...";
1291 SNESSetJacobian(snes,J,JPrec,SNESMF_form_jacobian_pn_2E,this);
1292 }
1293 else
1294 {
1295 SNESSetJacobian(snes,J,J,SNES_form_jacobian_pn_2E,this);
1296 }
1297
1298
1299 // set the newton method
1300 SNESSetType(snes,SNESLS); //default method
1301 // the maximum number of iterations
1302 PetscInt maxit = sv.maxit;
1303 if(sv.Type==EQUILIBRIUM) maxit = 10*sv.maxit;
1304
1305 if(sv.NS==LineSearch || sv.NS==Basic)
1306 {
1307 if(sv.NS==LineSearch)
1308 SNESLineSearchSet(snes,SNESLineSearchCubic,PETSC_NULL);
1309 if(sv.NS==Basic)
1310 SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);
1311
1312 if(sv.Damping==DampingBankRose)
1313 SNESLineSearchSetPostCheck(snes,LimitorByBankRose_L2E,this);
1314 else if(sv.Damping==DampingPotential)
1315 SNESLineSearchSetPostCheck(snes,LimitorByPotential_L2E,this);
1316 else
1317 SNESLineSearchSetPostCheck(snes,LimitorNonNegativeCarrier_L2E,this);
1318
1319 SNESSetTolerances(snes,1e-12*N,1e-14,1e-9,maxit,1000);
1320 SNESSetConvergenceTest(snes,LineSearch_ConvergenceTest_L2E,this);
1321 }
1322 else if(sv.NS==TrustRegion)
1323 {
1324 SNESSetType(snes,SNESTR);
1325 SNESSetTolerances(snes,1e-12*N,1e-14,1e-9,maxit,1000);
1326 SNESSetTrustRegionTolerance(snes,1e-20);
1327 SNESSetConvergenceTest(snes,TrustRegion_ConvergenceTest_L2E,this);
1328 //must set TR delta0 to a sufficient large value, or TR can't find real solution.
1329 SNES_TR *neP = (SNES_TR*)snes->data;
1330 neP->delta0 = N;
1331 }
1332
1333 //SNESMonitorSet(snes, SNESMonitorDefault,PETSC_NULL,PETSC_NULL);
1334
1335 SNESGetKSP(snes,&ksp);
1336 KSPGetPC(ksp,&pc);
1337 //KSPMonitorSet(ksp,KSPMonitorDefault,PETSC_NULL,PETSC_NULL);
1338 if(sv.LS=="lu"||sv.LS=="superlu"||sv.LS=="umfpack")
1339 {
1340 KSPSetType(ksp,KSPPREONLY);
1341 PCSetType(pc,PCLU);
1342 PCFactorSetReuseOrdering(pc,PETSC_TRUE);
1343 PCFactorSetPivoting(pc,1.0);
1344 PCFactorReorderForNonzeroDiagonal(pc,1e-14);
1345 PCFactorSetShiftNonzero(pc,1e-12);
1346 }
1347 else
1348 {
1349 KSPSetType(ksp,sv.LS.c_str());
1350 if(sv.LS=="gmres") KSPGMRESSetRestart(ksp,150);
1351 KSPSetTolerances(ksp,1e-12*N,1e-20*N,PETSC_DEFAULT,max(35,N/10));
1352 PCSetType(pc,PCILU);
1353 PCFactorSetLevels(pc,4);
1354 PCFactorSetShiftNonzero(pc,1e-14);
1355 PCFactorReorderForNonzeroDiagonal(pc,1e-14);
1356 }
1357 KSPSetFromOptions(ksp);
1358 SNESSetFromOptions(snes);
1359
1360 gss_log.string_buf()<<"done\n";
1361 gss_log.record();
1362 return 0;
1363 }
1364
1365
1366
1367 /* ----------------------------------------------------------------------------
1368 * DDM_Solver_L2E::do_solve: This function solve the problem
1369 */
1370 int DDM_Solver_L2E::do_solve(SolveDefine &sv)
1371 {
1372
1373 if(sv.Type==EQUILIBRIUM)
1374 solve_equ(sv);
1375
1376 else if(sv.Type==STEADYSTATE)
1377 solve_steadystate(sv);
1378
1379 else if(sv.Type==DCSWEEP)
1380 solve_dcsweep(sv);
1381
1382 else if(sv.Type==TRANSIENT)
1383 solve_transient(sv);
1384
1385 else if(sv.Type==TRACE)
1386 solve_iv_trace(sv);
1387
1388 else
1389 {
1390 gss_log.string_buf()<<"DDML2E not support this solver type!\n";
1391 gss_log.record();
1392 }
1393 return 0;
1394
1395 }
1396
1397
1398 /* ----------------------------------------------------------------------------
1399 * DDM_Solver_L2E::solve_equ: This function compute equilibrium
1400 * all the stimulate source(s) are zero. time step set to inf
1401 */
1402 int DDM_Solver_L2E::solve_equ(SolveDefine &sv)
1403 {
1404 gss_log.string_buf()<<"Compute equilibrium"<<"\n";
1405 gss_log.record();
1406 bc.Clear_Vapp();
1407 bc.Clear_Iapp();
1408 ODE_formula.TimeDependent = false;
1409 ODE_formula.dt = 1e100;
1410 ODE_formula.clock = 0.0;
1411 its = 0;
1412
1413 for(int i=0;i<zone_num;i++)
1414 if(zonedata[i]->material_type == Semiconductor)
1415 {
1416 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1417 pzonedata->HighFieldMobility= false;
1418 pzonedata->ImpactIonization = false;
1419 pzonedata->IIType = sv.IIType;
1420 pzonedata->BandBandTunneling = false;
1421 pzonedata->IncompleteIonization = false;
1422 pzonedata->QuantumMechanical = false;
1423 pzonedata->Fermi = sv.Fermi;
1424 pzonedata->EJModel = sv.EJModel;
1425 }
1426 // for equilibrium compute, low field mobility is ok.
1427 solver_pre_compute();
1428 probe_open(EQUILIBRIUM);
1429 SNESSolve(snes,PETSC_NULL,x);
1430 solution_update();
1431 probe(EQUILIBRIUM,0);
1432 probe_close();
1433 SNESGetConvergedReason(snes,&reason);
1434 if(reason<0)
1435 {
1436 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
1437 <<" "<<SNESConvergedReasons[reason]<<" !residual norm = "<<norm<<"\n\n\n";
1438 gss_log.record();
1439 }
1440 return 0;
1441 }
1442
1443
1444 /* ----------------------------------------------------------------------------
1445 * DDM_Solver_L2E::solve_steadystate: This function compute steadystate
1446 * set all the stimulate source(s) at there value of t=0. time step set to inf
1447 */
1448 int DDM_Solver_L2E::solve_steadystate(SolveDefine &sv)
1449 {
1450
1451 gss_log.string_buf()<<"Compute steady-state"<<"\n";
1452 gss_log.record();
1453 bc.Update_Vapp(0);
1454 bc.Update_Iapp(0);
1455 optgen_update(0);
1456 ODE_formula.TimeDependent = false;
1457 ODE_formula.dt = 1e100;
1458 its = 0;
1459
1460
1461 for(int i=0;i<zone_num;i++)
1462 if(zonedata[i]->material_type == Semiconductor)
1463 {
1464 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1465 pzonedata->HighFieldMobility= sv.HighFieldMobility;
1466 pzonedata->ImpactIonization = sv.ImpactIonization;
1467 pzonedata->IIType = sv.IIType;
1468 pzonedata->BandBandTunneling = sv.BandBandTunneling;
1469 pzonedata->IncompleteIonization = sv.IncompleteIonization;
1470 pzonedata->QuantumMechanical = sv.QuantumMechanical;
1471 pzonedata->Fermi = sv.Fermi;
1472 pzonedata->EJModel = sv.EJModel;
1473 }
1474 solver_pre_compute();
1475 probe_open(STEADYSTATE);
1476 SNESSolve(snes,PETSC_NULL,x);
1477 solution_update();
1478 probe(STEADYSTATE,0);
1479 probe_close();
1480 SNESGetConvergedReason(snes,&reason);
1481 if(reason<0)
1482 {
1483 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
1484 <<" "<<SNESConvergedReasons[reason]<<" !residual norm = "<<norm<<"\n\n\n";
1485 gss_log.record();
1486 }
1487 return 0;
1488
1489 }
1490
1491 /* ----------------------------------------------------------------------------
1492 * DDM_Solver_L2E::solve_dcsweep: This function compute IV curve
1493 * time step set to inf
1494 */
1495 int DDM_Solver_L2E::solve_dcsweep(SolveDefine &sv)
1496 {
1497 PetscScalar I=0;
1498 FILE *fiv;
1499 PetscScalar current_scale_mA = scale_unit.s_mA;
1500 PetscScalar voltage_scale_V =scale_unit.s_volt;
1501
1502 bc.Update_Vapp(0);
1503 bc.Update_Iapp(0);
1504 optgen_update(0);
1505 ODE_formula.TimeDependent = false;
1506 ODE_formula.dt = 1e100;
1507 ODE_formula.clock = 0.0;
1508
1509 //set which electrodes are required to record IV
1510 for(int i=0;i<zone_num;i++)
1511 {
1512 if(zonedata[i]->material_type == Semiconductor)
1513 {
1514 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1515 pzonedata->HighFieldMobility= sv.HighFieldMobility;
1516 pzonedata->ImpactIonization = sv.ImpactIonization;
1517 pzonedata->IIType = sv.IIType;
1518 pzonedata->BandBandTunneling = sv.BandBandTunneling;
1519 pzonedata->IncompleteIonization = sv.IncompleteIonization;
1520 pzonedata->QuantumMechanical = sv.QuantumMechanical;
1521 pzonedata->EJModel = sv.EJModel;
1522 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1523 for(int j=0;j<pzonedata->electrode.size();j++)
1524 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1525 {
1526 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1527 sv.Electrode_Record_Index.push_back(k);
1528 }
1529 }
1530 if(zonedata[i]->material_type == Insulator)
1531 {
1532 ISZone * pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1533 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1534 for(int j=0;j<pzonedata->electrode.size();j++)
1535 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1536 {
1537 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1538 sv.Electrode_Record_Index.push_back(k);
1539 }
1540 }
1541 }
1542
1543
1544 // output DC Scan information
1545 if(sv.Electrode_VScan!=-1)
1546 {
1547 if(!sv.Electrode_VScan_Name.size())
1548 {
1549 gss_log.string_buf()<<"No VScan Electrode Specified.\n";
1550 gss_log.record();
1551 return 1;
1552 }
1553 gss_log.string_buf()<<"DC voltage scan from "<<sv.VStart
1554 <<" step "<<sv.VStep
1555 <<" to "<<sv.VStop<<"\n";
1556 gss_log.record();
1557 }
1558 else
1559 {
1560 gss_log.string_buf()<<"DC current scan from "<<sv.IStart/current_scale_mA
1561 <<" step "<<sv.IStep/current_scale_mA
1562 <<" to "<<sv.IStop/current_scale_mA<<"\n";
1563 gss_log.record();
1564 }
1565
1566 // prepare IV file. If no file given, output to screen.
1567 if(!sv.IVFile.empty())
1568 {
1569 if(sv.IVFileAppend)
1570 {
1571 fiv=fopen(sv.IVFile.c_str(),"a");
1572 if(!fiv) fiv=stdout;
1573 }
1574 else
1575 {
1576 fiv=fopen(sv.IVFile.c_str(),"w");
1577 if(!fiv) fiv=stdout;
1578 fprintf(fiv,"#");
1579 for(int j=0;j<sv.Electrode_Record.size();j++)
1580 fprintf(fiv,"Vp(%s)[V] I(%s)[mA] ",
1581 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str(),
1582 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str()
1583 );
1584 fprintf(fiv,"\n");
1585 }
1586 }
1587 else
1588 fiv = stdout;
1589
1590 // begin
1591 if(sv.Electrode_VScan!=-1)
1592 {
1593 PetscScalar V = sv.VStart;
1594 PetscScalar VStep = sv.VStep;
1595 for(int i=0;i<sv.Electrode_VScan_Name.size();i++)
1596 bc.Set_electrode_type(sv.Electrode_VScan_Name[i].c_str(),VoltageBC);
1597 int diverged_retry=0;
1598 probe_open(DCSWEEP_VSCAN);
1599 do
1600 {
1601 its = 0;
1602 gss_log.string_buf()<<"DC Scan: V("<<sv.Electrode_VScan_Name[0];
1603 for(int i=1;i<sv.Electrode_VScan_Name.size();i++)
1604 gss_log.string_buf()<<", "<<sv.Electrode_VScan_Name[i];
1605 gss_log.string_buf()<<") = "<<V/voltage_scale_V<<" V"<<"\n";
1606 gss_log.record();
1607 for(int i=0;i<sv.Electrode_VScan_Name.size();i++)
1608 bc.Set_Vapp(sv.Electrode_VScan_Name[i].c_str(),V);
1609 solver_pre_compute();
1610 SNESSolve(snes,PETSC_NULL,x);
1611 SNESGetConvergedReason(snes,&reason);
1612
1613 if(reason>0 && norm==norm) //ok, converged.
1614 {
1615 diverged_retry=0;
1616 solution_update();
1617 probe(DCSWEEP_VSCAN,V);
1618 if(fabs(VStep) < fabs(sv.VStep))
1619 VStep *= 1.1;
1620 V+=VStep;
1621 if(V*sv.VStep > sv.VStop*sv.VStep && V*sv.VStep < (sv.VStop + VStep - 1e-10*VStep)*sv.VStep)
1622 V=sv.VStop;
1623 }
1624 else // oh, diverged... reduce step and try again
1625 {
1626 if(++diverged_retry>=8)
1627 {
1628 gss_log.string_buf()<<"------> Too many failed steps, give up tring.\n\n\n";
1629 gss_log.record();
1630 break;
1631 }
1632 diverged_recovery();
1633 VStep /= 2.0;
1634 V-=VStep;
1635 gss_log.string_buf()<<"------> nonlinear solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1636 gss_log.record();
1637 continue;
1638 }
1639
1640 for(int j=0;j<sv.Electrode_Record.size();j++)
1641 {
1642 I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new();
1643 fprintf(fiv,"%lf\t%e\t",
1644 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential_new()/voltage_scale_V),
1645 double(I/current_scale_mA));
1646 }
1647 if(sv.Electrode_Record.size())
1648 {
1649 fprintf(fiv,"\n");
1650 fflush(fiv);
1651 }
1652 }
1653 while(V*sv.VStep < (sv.VStop + 0.5*VStep)*sv.VStep);
1654 }
1655
1656 if(sv.Electrode_IScan!=-1)
1657 {
1658 bc.Set_electrode_type(sv.Electrode_IScan_Name.c_str(),CurrentBC);
1659 PetscScalar I = sv.IStart;
1660 PetscScalar IStep = sv.IStep;
1661 int diverged_retry=0;
1662 probe_open(DCSWEEP_ISCAN);
1663 do
1664 {
1665 its = 0;
1666 gss_log.string_buf()<<"DC Scan: I("<<sv.Electrode_IScan_Name<<") = "<<I/current_scale_mA<<" mA"<<"\n";
1667 gss_log.record();
1668 bc.Set_Iapp(sv.Electrode_IScan_Name.c_str(),I);
1669 solver_pre_compute();
1670 SNESSolve(snes,PETSC_NULL,x);
1671 SNESGetConvergedReason(snes,&reason);
1672
1673 if(reason>0 && norm==norm) //ok, converged.
1674 {
1675 diverged_retry=0;
1676 solution_update();
1677 probe(DCSWEEP_ISCAN,I);
1678 if(fabs(IStep) < fabs(sv.IStep))
1679 IStep *= 1.1;
1680 I+=IStep;
1681 if(I*sv.IStep > sv.IStop*sv.IStep && I*sv.IStep < (sv.IStop + IStep - 1e-10*IStep)*sv.IStep)
1682 I=sv.IStop;
1683 }
1684 else // oh, diverged... reduce step and try again
1685 {
1686 if(++diverged_retry>=8) //failed 8 times, stop tring
1687 {
1688 gss_log.string_buf()<<"------> Too many failed steps, give up tring.\n\n\n";
1689 gss_log.record();
1690 break;
1691 }
1692 diverged_recovery();
1693 IStep /= 2.0;
1694 I-=IStep;
1695 gss_log.string_buf()<<"------> nonlinear solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1696 gss_log.record();
1697 continue;
1698 }
1699
1700 for(int j=0;j<sv.Electrode_Record.size();j++)
1701 {
1702 fprintf(fiv,"%lf\t%e\t",
1703 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential_new()/voltage_scale_V),
1704 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new()/current_scale_mA));
1705 }
1706 if(sv.Electrode_Record.size())
1707 {
1708 fprintf(fiv,"\n");
1709 fflush(fiv);
1710 }
1711 }
1712 while(I*sv.IStep < (sv.IStop+0.5*IStep)*sv.IStep);
1713 }
1714
1715 if(!sv.IVFile.empty()) fclose(fiv);
1716 probe_close();
1717 return 0;
1718
1719 }
1720
1721
1722 /* ----------------------------------------------------------------------------
1723 * DDM_Solver_L2E::solve_iv_trace: This function use continuation method to trace
1724 * IV curve automatically
1725 */
1726 int DDM_Solver_L2E::solve_iv_trace(SolveDefine &sv)
1727 {
1728 int error=0;
1729 int first_step=1;
1730 int slope_flag=0;
1731 Vec pdI_pdw;
1732 Vec pdF_pdV;
1733 Vec pdw_pdV;
1734 PetscScalar pdI_pdV;
1735 PetscScalar dI_dV;
1736 PetscScalar I=0;
1737 PetscScalar V = sv.VStart;
1738 PetscScalar VStep = sv.VStep;
1739 PetscScalar r,Rload=0,Rload_new,Rref;
1740 PetscScalar angle,slope=PI/2,slope_new,slope_chord;
1741 FILE *fiv;
1742 PetscScalar current_scale_mA = scale_unit.s_mA;
1743 PetscScalar voltage_scale_V =scale_unit.s_volt;
1744 int bc_index;
1745 BaseBC * pbc = NULL;
1746 SMCZone *pz = NULL;
1747 int electrode_index;
1748 KSP kspc;
1749 PC pcc;
1750 VecCreateSeq(PETSC_COMM_SELF,N,&pdI_pdw);
1751 VecDuplicate(pdI_pdw,&pdF_pdV);
1752 VecDuplicate(pdI_pdw,&pdw_pdV);
1753 KSPCreate(MPI_COMM_SELF,&kspc);
1754 KSPSetType(kspc,KSPPREONLY);
1755 KSPGetPC(kspc,&pcc);
1756 PCSetType(pcc,PCLU);
1757 PCFactorReorderForNonzeroDiagonal(pcc,1e-14);
1758 PCFactorSetShiftNonzero(pcc,1e-12);
1759
1760 bc.Update_Vapp(0);
1761 bc.Update_Iapp(0);
1762 bc.Set_electrode_type(sv.Electrode_VScan_Name[0].c_str(),VoltageBC);
1763 bc.Set_Vapp(sv.Electrode_VScan_Name[0].c_str(),V);
1764 ODE_formula.TimeDependent = false;
1765 ODE_formula.dt = 1e100;
1766
1767 // output TRACE information
1768 gss_log.string_buf()<<"IV automatically trace by continuation method\n";
1769 gss_log.record();
1770
1771 //set which electrodes are required to do trace or record IV
1772 for(int i=0;i<zone_num;i++)
1773 {
1774 if(zonedata[i]->material_type == Semiconductor)
1775 {
1776 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1777 pzonedata->HighFieldMobility= sv.HighFieldMobility;
1778 pzonedata->ImpactIonization = sv.ImpactIonization;
1779 pzonedata->IIType = sv.IIType;
1780 pzonedata->BandBandTunneling = sv.BandBandTunneling;
1781 pzonedata->IncompleteIonization = sv.IncompleteIonization;
1782 pzonedata->QuantumMechanical = sv.QuantumMechanical;
1783 pzonedata->Fermi = sv.Fermi;
1784 pzonedata->EJModel = sv.EJModel;
1785 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1786 for(int j=0;j<pzonedata->electrode.size();j++)
1787 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1788 {
1789 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1790 sv.Electrode_Record_Index.push_back(k);
1791 }
1792 for(int j=0;j<pzonedata->electrode.size();j++)
1793 {
1794 if(bc.is_electrode_label(pzonedata->electrode[j],sv.Electrode_VScan_Name[0].c_str()))
1795 {
1796 pz = pzonedata;
1797 electrode_index = j;
1798 bc_index = pzonedata->electrode[j];
1799 pbc = bc.Get_pointer(bc_index);
1800 }
1801 }
1802 }
1803 if(zonedata[i]->material_type == Insulator)
1804 {
1805 ISZone * pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1806 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1807 for(int j=0;j<pzonedata->electrode.size();j++)
1808 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1809 {
1810 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1811 sv.Electrode_Record_Index.push_back(k);
1812 }
1813 }
1814 }
1815
1816 //check error
1817 if(!pbc || !pz )
1818 {
1819 gss_log.string_buf()<<"I can't find suitable electrode to do IV trace.\n";
1820 gss_log.record();
1821 error = 1;
1822 goto trace_end;
1823 }
1824 if(pbc->BCType!=OhmicContact && pbc->BCType!=SchottkyContact)
1825 {
1826 gss_log.string_buf()<<"IV trace can only be performanced to Ohmic or Schottky contact.\n";
1827 gss_log.record();
1828 error = 1;
1829 goto trace_end;
1830 }
1831
1832 // prepare IV file. If no file given, output to screen.
1833 if(!sv.IVFile.empty())
1834 {
1835 fiv=fopen(sv.IVFile.c_str(),"w");
1836 fprintf(fiv,"#");
1837 for(int j=0;j<sv.Electrode_Record.size();j++)
1838 fprintf(fiv,"Vp(%s)[V] I(%s)[mA] ",
1839 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str(),
1840 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str()
1841 );
1842 fprintf(fiv,"\n");
1843 }
1844 else
1845 fiv = stdout;
1846
1847 // initial condition check
1848 pbc->Set_R(Rload);
1849 solver_pre_compute();
1850 SNESSolve(snes,PETSC_NULL,x);
1851 SNESGetConvergedReason(snes,&reason);
1852 if(reason<0)
1853 {
1854 gss_log.string_buf()<<"I can't get convergence even at initial point, please give a good initial condition.\n\n";
1855 gss_log.record();
1856 error = 1;
1857 goto trace_end;
1858 }
1859 solution_update();
1860 // output IV curve
1861 for(int j=0;j<sv.Electrode_Record.size();j++)
1862 {
1863 I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new();
1864 fprintf(fiv,"%lf\t%e\t",
1865 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential()/voltage_scale_V),
1866 double(I/current_scale_mA));
1867 }
1868 if(sv.Electrode_Record.size())
1869 {
1870 fprintf(fiv,"\n");
1871 fflush(fiv);
1872 }
1873
1874 //loop here
1875 while(pbc->Get_Potential()*sv.VStep < sv.VStop*sv.VStep && fabs(pbc->Get_Current())<sv.IStop)
1876 {
1877 // for last Rload, increase VStep
1878 gss_log.string_buf()<<"Trace for V("<<sv.Electrode_VScan_Name[0]<<")="<<pbc->Get_Potential()<<"(V)\n";
1879 gss_log.record();
1880
1881 int recovery=0;
1882 do
1883 {
1884 if(!slope_flag && first_step)
1885 V+=VStep;
1886 if(!slope_flag && !first_step)
1887 V*=1+VStep*cos(slope)/pbc->Get_Potential();
1888 bc.Set_Vapp(sv.Electrode_VScan_Name[0].c_str(),V);
1889 SNESSolve(snes,PETSC_NULL,x);
1890 SNESGetConvergedReason(snes,&reason);
1891 if(reason<0)
1892 {
1893 gss_log.string_buf()<<"I can't get convergence at this step, do recovery...\n\n";
1894 gss_log.record();
1895 diverged_recovery();
1896 V/=1+VStep*cos(slope)/pbc->Get_Potential();
1897 VStep/=2;
1898 recovery++;
1899 }
1900 if(recovery>8)
1901 {
1902 gss_log.string_buf()<<"Too many failure step, give up tring.\n\n\n";
1903 gss_log.record();
1904 error = 1;
1905 goto trace_end;
1906 }
1907 }
1908 while(reason<0);
1909
1910 // for the new bias point, recompute slope
1911 // calculate the dynamic resistance of IV curve by different approximation
1912 //r = (pbc->Get_Potential_new()-pbc->Get_Potential())/(pbc->Get_Current_new()-pbc->Get_Current());
1913
1914 if(pbc->BCType==OhmicContact)
1915 pz->F2E_om_electrode_Trace(electrode_index,pdI_pdV,pdI_pdw,pdF_pdV,x,&JTmp,&J,ODE_formula,zofs,bc,DeviceDepth);
1916 if(pbc->BCType==SchottkyContact)
1917 pz->F2E_stk_electrode_Trace(electrode_index,pdI_pdV,pdI_pdw,pdF_pdV,x,&JTmp,&J,ODE_formula,zofs,bc,DeviceDepth);
1918 KSPSetOperators(kspc,J,J,SAME_NONZERO_PATTERN);
1919 KSPSetUp(kspc);
1920 KSPSolve(kspc,pdF_pdV,pdw_pdV);
1921 VecDot(pdI_pdw,pdw_pdV,&dI_dV);
1922 r=1/dI_dV;
1923 //printf("dI/dV=%e pdI/pdV=%e VStep=%e\n",1/r,dI_dV,VStep);
1924
1925 if(first_step)
1926 {
1927 Rref = pbc->Get_Potential_new()/pbc->Get_Current_new();
1928 slope=slope_new=atan(Rref/r);
1929 }
1930 else
1931 {
1932 Rref = pbc->Get_Potential()/pbc->Get_Current(); // Rref, for scaling the dynamic resistance
1933 slope_new = atan(Rref/r);
1934 }
1935
1936 //resolve turning point
1937 if(tan(slope)*tan(slope_new)<0) //the sign is changed
1938 {
1939 PetscScalar chord = (pbc->Get_Potential_new()-pbc->Get_Potential())/(pbc->Get_Current_new()-pbc->Get_Current());
1940 slope_chord = atan(Rref/chord);
1941 //printf("slope_chord=%e slope=%e slope_new=%e \n",slope_chord/PI*180, slope/PI*180, slope_new/PI*180);
1942 //for vertical turning point
1943 if(fabs(slope)/PI*180 >70 && fabs(slope_new)/PI*180 >70 &&
1944 ((slope_chord*slope>0 && fabs(slope_chord)>fabs(slope)) || (slope_chord*slope_new>0 && fabs(slope_chord)>fabs(slope_new))))
1945 {
1946 VStep*=-1.0;
1947 //gss_log.string_buf()<<"vertical turning point...\n";
1948 //gss_log.record();
1949 }
1950 //for horizontal turning point
1951 if( fabs(slope)/PI*180 <20 && fabs(slope_new)/PI*180 <20 &&
1952 ((slope_chord*slope>0 && fabs(slope_chord)<fabs(slope)) || (slope_chord*slope_new>0 && fabs(slope_chord)<fabs(slope_new))))
1953 {
1954 VStep*=1.0;
1955 //gss_log.string_buf()<<"horizontal turning point...\n";
1956 //gss_log.record();
1957 }
1958 }
1959
1960 // check the slope change
1961 angle=fabs(slope-slope_new);
1962 if(angle>PI/2) angle=PI-angle;
1963 if(angle<PI/36) VStep*=1.5; // slope change less than 5 degree
1964 else if(angle<PI/24) VStep=VStep; // slope change less than 10 degree, but greater than 5 degree
1965 else if(angle<PI/12) VStep/=2; // slope change less than 15 degree, but greater than 10 degree
1966 else // slope greater than 15 degree, reject this solution
1967 {
1968 //printf("slope change old=%e new=%e angle=%e \n",slope/PI*180, slope_new/PI*180, angle/PI*180);
1969 //printf("V=%e VStep=%e\n",V,VStep);
1970 gss_log.string_buf()<<"Slope of IV curve changes too quickly, do recovery...\n\n";
1971 gss_log.record();
1972 V/=1+VStep*cos(slope)/pbc->Get_Potential();
1973 VStep/=2;
1974 slope_flag=1;
1975 //printf("V=%e VStep=%e\n",V,VStep);
1976 diverged_recovery();
1977 continue;
1978 }
1979 if(fabs(VStep*cos(slope))>sv.VStepMax) VStep=VStep/fabs(VStep)*sv.VStepMax/cos(slope);
1980 // ok, update solutions
1981 solution_update();
1982 slope_flag=0;
1983 // since Rload will be changed, we should change bias voltage to meet the truncation point.
1984
1985 Rload_new=Rref/r;
1986 V+=pbc->Get_Current()*(Rload_new-Rload);
1987 pbc->Set_R(Rload_new);
1988 bc.Set_Vapp(sv.Electrode_VScan_Name[0].c_str(),V);
1989 Rload=Rload_new;
1990 slope = slope_new;
1991 // output IV curve
1992 for(int j=0;j<sv.Electrode_Record.size();j++)
1993 {
1994 I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new();
1995 fprintf(fiv,"%lf\t%e\t",
1996 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential()/voltage_scale_V),
1997 double(I/current_scale_mA));
1998 }
1999 if(sv.Electrode_Record.size())
2000 {
2001 fprintf(fiv,"\n");
2002 fflush(fiv);
2003 }
2004
2005 first_step=0;
2006 }
2007
2008 trace_end:
2009 KSPDestroy(kspc);
2010 VecDestroy(pdI_pdw);
2011 VecDestroy(pdF_pdV);
2012 VecDestroy(pdw_pdV);
2013 if(!sv.IVFile.empty()) fclose(fiv);
2014 return error;
2015 }
2016
2017
2018 void DDM_Solver_L2E::LET_norm_estimat(PetscScalar & r)
2019 {
2020 PetscScalar hn= ODE_formula.dt;
2021 PetscScalar hn1= ODE_formula.dt_last;
2022 PetscScalar hn2= ODE_formula.dt_last_last;
2023 PetscScalar eps_r=1e-3,eps_a=1e-4;
2024
2025 VecZeroEntries(xp);
2026 VecZeroEntries(LTE);
2027
2028 if(ODE_formula.BDF_Type==BDF1)
2029 {
2030 VecAXPY(xp,1+hn/hn1,x_n);
2031 VecAXPY(xp,-hn/hn1,x_n1);
2032 VecAXPY(LTE,hn/(hn+hn1),x);
2033 VecAXPY(LTE,-hn/(hn+hn1),xp);
2034 }
2035 else if(ODE_formula.BDF_Type==BDF2)
2036 {
2037 PetscScalar cn=1+hn*(hn+2*hn1+hn2)/(hn1*(hn1+hn2));
2038 PetscScalar cn1=-hn*(hn+hn1+hn2)/(hn1*hn2);
2039 PetscScalar cn2=hn*(hn+hn1)/(hn2*(hn1+hn2));
2040
2041 VecAXPY(xp,cn,x_n);
2042 VecAXPY(xp,cn1,x_n1);
2043 VecAXPY(xp,cn2,x_n2);
2044 VecAXPY(LTE, hn/(hn+hn1+hn2),x);
2045 VecAXPY(LTE,-hn/(hn+hn1+hn2),xp);
2046 }
2047
2048 PetscScalar *xx,*ll;
2049 VecGetArray(x,&xx);
2050 VecGetArray(LTE,&ll);
2051 int offset=0;
2052 int N=0;
2053 for(int z=0;z<zone_num;z++)
2054 {
2055 if(zonedata[z]->material_type==Semiconductor)
2056 {
2057 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
2058 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
2059 {
2060 ll[offset+0]=0;
2061 ll[offset+1]=ll[offset+1]/(eps_r*xx[offset+1]+eps_a);
2062 ll[offset+2]=ll[offset+2]/(eps_r*xx[offset+2]+eps_a);
2063 ll[offset+3]=ll[offset+3]/(eps_r*xx[offset+3]+eps_a);
2064 offset += 4;
2065 N+=3;
2066 }
2067 for(int i=0;i<pzonedata->electrode.size();i++)
2068 {
2069 ll[offset+0]=0;
2070 offset += 1;
2071 }
2072 }
2073 else if(zonedata[z]->material_type==Insulator)
2074 {
2075 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
2076 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
2077 {
2078 ll[offset+0]=0;
2079 ll[offset+1]=ll[offset+1]/(eps_r*xx[offset+1]+eps_a);
2080 offset += 2;
2081 N+=1;
2082 }
2083 for(int i=0;i<pzonedata->electrode.size();i++)
2084 {
2085 ll[offset+0]=0;
2086 offset += 1;
2087 }
2088 }
2089 else if(zonedata[z]->material_type==Conductor)
2090 {
2091 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
2092 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
2093 {
2094 ll[offset+0]=0;
2095 ll[offset+1]=ll[offset+1]/(eps_r*xx[offset+1]+eps_a);
2096 offset += 2;
2097 N+=1;
2098 }
2099 }
2100 }
2101 VecRestoreArray(x,&xx);
2102 VecRestoreArray(LTE,&ll);
2103 VecNorm(LTE,NORM_2,&r);
2104
2105 r/=sqrt(N);
2106 }
2107
2108
2109 /* ----------------------------------------------------------------------------
2110 * DDM_Solver_L2E::solve_transient: This function does transient simulation.
2111 */
2112 int DDM_Solver_L2E::solve_transient(SolveDefine &sv)
2113 {
2114
2115 PetscScalar I=0;
2116 FILE *fiv;
2117
2118 PetscScalar current_scale_mA = scale_unit.s_mA;
2119 PetscScalar voltage_scale_V =scale_unit.s_volt;
2120 int diverged_retry=0;
2121 int step=0;
2122 PetscScalar r=1.0;
2123 clock = sv.TStart+sv.TStep;
2124 ODE_formula.TimeDependent = true;
2125 ODE_formula.BDF_Type = sv.BDF_Type;
2126 ODE_formula.dt = sv.TStep;
2127 if(sv.BDF_Type==BDF2)
2128 {
2129 ODE_formula.BDF2_restart = true;
2130 }
2131
2132 for(int i=0;i<zone_num;i++)
2133 {
2134 if(zonedata[i]->material_type == Semiconductor)
2135 {
2136 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
2137 pzonedata->HighFieldMobility= sv.HighFieldMobility;
2138 pzonedata->ImpactIonization = sv.ImpactIonization;
2139 pzonedata->IIType = sv.IIType;
2140 pzonedata->BandBandTunneling = sv.BandBandTunneling;
2141 pzonedata->IncompleteIonization = sv.IncompleteIonization;
2142 pzonedata->QuantumMechanical = sv.QuantumMechanical;
2143 pzonedata->Fermi = sv.Fermi;
2144 pzonedata->EJModel = sv.EJModel;
2145 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
2146 for(int j=0;j<pzonedata->electrode.size();j++)
2147 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
2148 {
2149 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
2150 sv.Electrode_Record_Index.push_back(k);
2151 }
2152 }
2153 if(zonedata[i]->material_type == Insulator)
2154 {
2155 ISZone * pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
2156 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
2157 for(int j=0;j<pzonedata->electrode.size();j++)
2158 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
2159 {
2160 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
2161 sv.Electrode_Record_Index.push_back(k);
2162 }
2163 }
2164 }
2165
2166 gss_log.string_buf()<<"Transient compute from "<<sv.TStart<<" ps step "<<sv.TStep<<" ps to "<<sv.TStop<<" ps"<<"\n";
2167 gss_log.record();
2168
2169 //open record file
2170 if(!sv.IVFile.empty())
2171 {
2172 if(sv.IVFileAppend)
2173 {
2174 fiv=fopen(sv.IVFile.c_str(),"a");
2175 if(!fiv) fiv=stdout;
2176 }
2177 else
2178 {
2179 fiv=fopen(sv.IVFile.c_str(),"w");
2180 if(!fiv) fiv=stdout;
2181 fprintf(fiv,"#");
2182 fprintf(fiv,"time[ps] ");
2183 for(int j=0;j<sv.Electrode_Record.size();j++)
2184 fprintf(fiv,"Vp(%s)[V] I(%s)[mA] ",
2185 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str(),
2186 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str()
2187 );
2188 fprintf(fiv,"\n");
2189 }
2190 }
2191 else
2192 fiv = stdout;
2193
2194 probe_open(TRANSIENT);
2195
2196 // the main loop of transient solver.
2197 do
2198 {
2199 gss_log.string_buf()<<"t = "<<clock<<" ps"<<"\n";
2200 gss_log.record();
2201 ODE_formula.clock = clock;
2202
2203 its = 0; //clear nonlinear iteration counter
2204 //update sources
2205 bc.Update_Vapp(clock);
2206 bc.Update_Iapp(clock);
2207 optgen_update(clock);
2208
2209 //we do solve here!
2210 solver_pre_compute();
2211 SNESSolve(snes,PETSC_NULL,x);
2212 SNESGetConvergedReason(snes,&reason);
2213
2214 //nonlinear solution diverged?
2215 if(reason<0 || !(norm==norm))
2216 {
2217 if(++diverged_retry>=8) //failed 8 times, stop tring
2218 {
2219 gss_log.string_buf()<<"------> Too many failed steps, give up tring.\n\n\n";
2220 gss_log.record();
2221 break;
2222 }
2223 diverged_recovery();
2224 gss_log.string_buf()<<"------> nonlinear solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
2225 gss_log.record();
2226 ODE_formula.dt/=2.0; // reduce time step by a factor of two
2227 clock-=ODE_formula.dt;
2228 if(clock<sv.TStart) clock=sv.TStart;
2229 continue;
2230 }
2231
2232 //ok, nonlinear solution converged.
2233 diverged_retry=0;
2234
2235 //do LTE estimation and auto time step control
2236 if(sv.AutoStep && ((ODE_formula.BDF_Type==BDF1 && step>=3) || (ODE_formula.BDF_Type==BDF2 && step>=4)))
2237 {
2238 LET_norm_estimat(r);
2239
2240 if(ODE_formula.BDF_Type==BDF1)
2241 r = pow(r,PetscScalar(-1.0/2));
2242 else if(ODE_formula.BDF_Type==BDF2)
2243 r = pow(r,PetscScalar(-1.0/3));
2244
2245 if(r<0.9) //reject this solution
2246 {
2247 diverged_recovery();
2248 gss_log.string_buf()<<"------> LTE too large, time step rejected...\n\n\n";
2249 gss_log.record();
2250 clock-=ODE_formula.dt;
2251 ODE_formula.dt*=r; // reduce time step by a factor of r
2252 clock+=ODE_formula.dt;
2253 continue;
2254 }
2255 ODE_formula.dt=ODE_formula.dt * (r > 2.0 ? 2.0 : r);
2256 if(ODE_formula.dt>4*sv.TStep) ODE_formula.dt=4*sv.TStep;
2257 }
2258 else
2259 {
2260 if(fabs(ODE_formula.dt) < fabs(sv.TStep))
2261 ODE_formula.dt*=1.1;
2262 }
2263
2264 //update solution
2265 step++;
2266 solution_update();
2267 VecCopy(x_n1,x_n2);
2268 VecCopy(x_n,x_n1);
2269 VecCopy(x,x_n);
2270 ODE_formula.dt_last_last = ODE_formula.dt_last;
2271 ODE_formula.dt_last = ODE_formula.dt;
2272
2273 //predict next solution
2274 if(sv.Predict && (ODE_formula.BDF_Type==BDF1 && step>=3) || (ODE_formula.BDF_Type==BDF2 && step>=4))
2275 {
2276 PetscScalar hn= ODE_formula.dt;
2277 PetscScalar hn1= ODE_formula.dt_last;
2278 PetscScalar hn2= ODE_formula.dt_last_last;
2279 VecZeroEntries(x);
2280 if(ODE_formula.BDF_Type==BDF1)
2281 {
2282 VecAXPY(x,1+hn/hn1,x_n);
2283 VecAXPY(x,-hn/hn1,x_n1);
2284 }
2285 else if(ODE_formula.BDF_Type==BDF2)
2286 {
2287 PetscScalar cn=1+hn*(hn+2*hn1+hn2)/(hn1*(hn1+hn2));
2288 PetscScalar cn1=-hn*(hn+hn1+hn2)/(hn1*hn2);
2289 PetscScalar cn2=hn*(hn+hn1)/(hn2*(hn1+hn2));
2290
2291 VecAXPY(x,cn,x_n);
2292 VecAXPY(x,cn1,x_n1);
2293 VecAXPY(x,cn2,x_n2);
2294 }
2295 }
2296
2297 //save solution to file
2298 probe(TRANSIENT,clock);
2299 if(sv.Electrode_Record.size())
2300 {
2301 fprintf(fiv,"%e\t",double(clock));
2302 for(int j=0;j<sv.Electrode_Record.size();j++)
2303 {
2304 I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new();
2305 fprintf(fiv,"%lf\t%e\t",
2306 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential_new()/voltage_scale_V),
2307 double(I/current_scale_mA));
2308 }
2309 fprintf(fiv,"%lf\t",double(r));
2310 if(sv.Electrode_Record.size())
2311 {
2312 fprintf(fiv,"\n");
2313 fflush(fiv);
2314 }
2315 }
2316
2317 //make sure we can terminat near TStop, the relative error should less than 1e-10.
2318 clock+=ODE_formula.dt;
2319 if(clock > sv.TStop && clock < (sv.TStop + ODE_formula.dt - 1e-10*ODE_formula.dt))
2320 {
2321 ODE_formula.dt -= clock - sv.TStop;
2322 clock=sv.TStop;
2323 }
2324
2325 //clear the first step flag of BDF2
2326 if(ODE_formula.BDF_Type==BDF2)
2327 ODE_formula.BDF2_restart = false;
2328
2329
2330 }while(clock<sv.TStop+0.5*ODE_formula.dt);
2331
2332 //close record file
2333 if(!sv.IVFile.empty())
2334 fclose(fiv);
2335 probe_close();
2336
2337 return 0;
2338
2339 }
2340
2341
2342
2343 /* ----------------------------------------------------------------------------
2344 * DDM_Solver_L2E::solver_pre_compute: This function precompute some parameters before
2345 * each newton iteration.
2346 */
2347 void DDM_Solver_L2E::solver_pre_compute()
2348 {
2349 for(int z=0;z<zone_num;z++)
2350 {
2351 if(zonedata[z]->material_type==Semiconductor)
2352 {
2353 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
2354 for(int i=0;i<zone[z].davcell.size();i++)
2355 {
2356 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
2357 SemiData *pfs = &(pzonedata->fs[i]);
2358 SemiAuxData *paux = &(pzonedata->aux[i]);
2359 pzonedata->mt->mapping(&zone[z].danode[i],paux,0);
2360 pzonedata->aux[i].mun = pzonedata->mt->mob->ElecMob(pfs->p,pfs->n,pfs->T,0,0,pfs->T);
2361 pzonedata->aux[i].mup = pzonedata->mt->mob->HoleMob(pfs->p,pfs->n,pfs->T,0,0,pfs->T);
2362 pzonedata->aux[i].affinity = pzonedata->mt->basic->Affinity(pfs->T);
2363 pzonedata->aux[i].density = pzonedata->mt->basic->Density(pfs->T);
2364 pzonedata->aux[i].Eg = pzonedata->mt->band->Eg(pfs->T);
2365 pzonedata->aux[i].Nc = pzonedata->mt->band->Nc(pfs->T);
2366 pzonedata->aux[i].Nv = pzonedata->mt->band->Nv(pfs->T);
2367 }
2368 }
2369 }
2370 }
2371
2372
2373
2374 /* ----------------------------------------------------------------------------
2375 * DDM_Solver_L2E::solution_update: This function restore solution data from SNES
2376 */
2377 void DDM_Solver_L2E::solution_update()
2378 {
2379 //------------------------------------------------------------
2380 PetscScalar *xx;
2381 VecGetArray(x,&xx);
2382 int offset=0;
2383 for(int z=0;z<zone_num;z++)
2384 {
2385 if(zonedata[z]->material_type==Semiconductor)
2386 {
2387 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
2388 PetscScalar kb = pzonedata->mt->kb;
2389 PetscScalar e = pzonedata->mt->e;
2390 for(int i=0;i<zone[z].davcell.size();i++)
2391 {
2392
2393 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
2394 pzonedata->fs[i].n_last = pzonedata->fs[i].n;
2395 pzonedata->fs[i].p_last = pzonedata->fs[i].p;
2396 pzonedata->fs[i].T_last = pzonedata->fs[i].T;
2397 pzonedata->fs[i].P = xx[offset+0];
2398 pzonedata->fs[i].n = fabs(xx[offset+1]);
2399 pzonedata->fs[i].p = fabs(xx[offset+2]);
2400 pzonedata->fs[i].T = xx[offset+3];
2401
2402 pzonedata->aux[i].affinity = pzonedata->mt->basic->Affinity(xx[offset+3]);
2403 pzonedata->aux[i].Eg = pzonedata->mt->band->Eg(xx[offset+3]);
2404 pzonedata->aux[i].Nc = pzonedata->mt->band->Nc(xx[offset+3]);
2405 pzonedata->aux[i].Nv = pzonedata->mt->band->Nv(xx[offset+3]);
2406 pzonedata->mt->mapping(&pzonedata->pzone->danode[i],&pzonedata->aux[i],0);
2407 PetscScalar nie = pzonedata->mt->band->nie(pzonedata->fs[i].T);
2408 pzonedata->aux[i].phi_intrinsic = pzonedata->fs[i].P + pzonedata->aux[i].affinity +
2409 kb*pzonedata->fs[i].T/e*log(pzonedata->aux[i].Nc/nie);
2410 pzonedata->aux[i].phin = pzonedata->aux[i].phi_intrinsic - log(fabs(pzonedata->fs[i].n)/nie)*kb*pzonedata->fs[i].T/e;
2411 pzonedata->aux[i].phip = pzonedata->aux[i].phi_intrinsic + log(fabs(pzonedata->fs[i].p)/nie)*kb*pzonedata->fs[i].T/e;
2412
2413 offset += 4;
2414 }
2415 pzonedata->F2E_efield_update(xx,zofs,bc,zonedata);
2416 for(int i=0;i<pzonedata->electrode.size();i++)
2417 {
2418 PetscScalar P = xx[offset];
2419 PetscScalar Pn = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
2420 PetscScalar C = bc.Get_pointer(pzonedata->electrode[i])->Get_C();
2421 PetscScalar I = bc.Get_pointer(pzonedata->electrode[i])->Get_Current_new();
2422 bc.Get_pointer(pzonedata->electrode[i])->Set_Cap_Current(C*(P-Pn)/ODE_formula.dt);
2423 bc.Get_pointer(pzonedata->electrode[i])->Set_Current(I);
2424 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential(P);
2425 offset += 1;
2426 }
2427 }
2428 if(zonedata[z]->material_type==Insulator)
2429 {
2430 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
2431 pzonedata->F2_efield_update(xx,zofs,bc,zonedata);
2432 for(int i=0;i<zone[z].davcell.size();i++)
2433 {
2434 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
2435 pzonedata->fs[i].T_last = pzonedata->fs[i].T;
2436 pzonedata->fs[i].P = xx[offset+0];
2437 pzonedata->fs[i].T = xx[offset+1];
2438 offset += 2;
2439 }
2440 for(int i=0;i<pzonedata->electrode.size();i++)
2441 {
2442 PetscScalar P = xx[offset];
2443 PetscScalar Pn = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
2444 PetscScalar C = bc.Get_pointer(pzonedata->electrode[i])->Get_C();
2445 PetscScalar I = bc.Get_pointer(pzonedata->electrode[i])->Get_Current_new();
2446 bc.Get_pointer(pzonedata->electrode[i])->Set_Cap_Current(C*(P-Pn)/ODE_formula.dt);
2447 bc.Get_pointer(pzonedata->electrode[i])->Set_Current(I);
2448 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential(P);
2449 offset += 1;
2450 }
2451 }
2452 if(zonedata[z]->material_type==Conductor)
2453 {
2454 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
2455 for(int i=0;i<zone[z].davcell.size();i++)
2456 {
2457 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
2458 pzonedata->fs[i].T_last = pzonedata->fs[i].T;
2459 pzonedata->fs[i].P = xx[offset+0];
2460 pzonedata->fs[i].T = xx[offset+1];
2461 offset += 2;
2462 }
2463 }
2464 }
2465 ODE_formula.dt_last = ODE_formula.dt;
2466
2467 VecRestoreArray(x,&xx);
2468 }
2469
2470 /* ----------------------------------------------------------------------------
2471 * DDM_Solver_L2E::diverged_recovery: This function recovery latest solution data
2472 * if SNES diverged.
2473 */
2474 void DDM_Solver_L2E::diverged_recovery()
2475 {
2476 //------------------------------------------------------------
2477 PetscScalar *xx;
2478 VecGetArray(x,&xx);
2479 int offset=0;
2480 for(int z=0;z<zone_num;z++)
2481 {
2482 if(zonedata[z]->material_type==Semiconductor)
2483 {
2484 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
2485 for(int i=0;i<zone[z].davcell.size();i++)
2486 {
2487 xx[offset+0] = pzonedata->fs[i].P;
2488 xx[offset+1] = pzonedata->fs[i].n;
2489 xx[offset+2] = pzonedata->fs[i].p;
2490 xx[offset+3] = pzonedata->fs[i].T;
2491 offset += 4;
2492 }
2493 for(int i=0;i<pzonedata->electrode.size();i++)
2494 {
2495 PetscScalar P = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
2496 xx[offset] = P;
2497 offset += 1;
2498 }
2499 }
2500 if(zonedata[z]->material_type==Insulator)
2501 {
2502 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
2503 for(int i=0;i<zone[z].davcell.size();i++)
2504 {
2505 xx[offset+0] = pzonedata->fs[i].P;
2506 xx[offset+1] = pzonedata->fs[i].T;
2507 offset += 2;
2508 }
2509 for(int i=0;i<pzonedata->electrode.size();i++)
2510 {
2511 PetscScalar P = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
2512 xx[offset] = P;
2513 offset += 1;
2514 }
2515 }
2516 if(zonedata[z]->material_type==Conductor)
2517 {
2518 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
2519 for(int i=0;i<zone[z].davcell.size();i++)
2520 {
2521 xx[offset+0] = pzonedata->fs[i].P;
2522 xx[offset+1] = pzonedata->fs[i].T;
2523 offset += 2;
2524 }
2525 }
2526 }
2527 VecRestoreArray(x,&xx);
2528 }
2529
2530 /* ----------------------------------------------------------------------------
2531 * DDM_Solver_L2E::destroy_solver: This function do destroy the nonlinear solver
2532 */
2533 int DDM_Solver_L2E::destroy_solver(SolveDefine &sv)
2534 {
2535 // free work space
2536 N = 0;
2537 zofs.clear();
2538 VecDestroy(x);
2539 VecDestroy(r);
2540 VecDestroy(x_n);
2541 VecDestroy(x_n1);
2542 VecDestroy(x_n2);
2543 VecDestroy(xp);
2544 VecDestroy(LTE);
2545 MatDestroy(J);
2546 MatDestroy(JTmp);
2547 if(mf_flg)
2548 MatDestroy(JPrec);
2549 SNESDestroy(snes);
2550 return 0;
2551 }
2552
General-purpose Semiconductor Simulator