?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fix portability.
[gss-tcad.git] / src / solver / mix1 / mix1.cc
blob8b89ac136d3d6efb47e6114e32d95c378e1a9717
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 29, 2007 */
14 /* */
15 /* Gong Ding */
16 /* gdiso@ustc.edu */
17 /* NINT, No.69 P.O.Box, Xi'an City, China */
18 /* */
19 /*****************************************************************************/
20 /* the interface data structure with ng-spice */
21
22 #include "mathfunc.h"
23 #include "mix1.h"
24 #include "private/kspimpl.h"
25 #include "private/snesimpl.h"
26 #include "src/snes/impls/tr/tr.h"
27 #include "log.h"
28 #ifdef HAVE_UNISTD_H
29 #include <unistd.h>
30 #endif
31
32 /* ----------------------------------------------------------------------------
33 * error_norm_pn_Mix1: This function compute X and RHS error norms.
34 */
35 void DDM_Mix_Solver_L1E:: error_norm_pn_Mix1(PetscScalar *x,PetscScalar *f)
36 {
37 // do clear
38 potential_norm = 0;
39 electron_norm = 0;
40 hole_norm = 0;
41
42 possion_norm = 0;
43 elec_continuty_norm = 0;
44 hole_continuty_norm = 0;
45
46 int offset=0;
47 for(int z=0;z<zone_num;z++)
48 {
49 if(zonedata[z]->material_type==Semiconductor)
50 {
51 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
52 for(int i=0;i<zone[z].davcell.size();i++)
53 {
54 potential_norm += x[offset+0]*x[offset+0];
55 electron_norm += x[offset+1]*x[offset+1];
56 hole_norm += x[offset+2]*x[offset+2];
57 possion_norm += f[offset+0]*f[offset+0];
58 elec_continuty_norm += f[offset+1]*f[offset+1];
59 hole_continuty_norm += f[offset+2]*f[offset+2];
60 offset += 3;
61 }
62 }
63 if(zonedata[z]->material_type==Insulator)
64 {
65 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
66 for(int i=0;i<zone[z].davcell.size();i++)
67 {
68 potential_norm += x[offset]*x[offset];
69 possion_norm += f[offset]*f[offset];
70 offset += 1;
71 }
72 for(int i=0;i<pzonedata->electrode.size();i++)
73 {
74 potential_norm += x[offset]*x[offset];
75 offset += 1;
76 }
77 }
78 if(zonedata[z]->material_type==Conductor)
79 {
80 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
81 for(int i=0;i<zone[z].davcell.size();i++)
82 {
83 potential_norm += x[offset]*x[offset];
84 possion_norm += f[offset]*f[offset];
85 offset += 1;
86 }
87 }
88 }
89
90 potential_norm = sqrt(potential_norm);
91 electron_norm = sqrt(electron_norm);
92 hole_norm = sqrt(hole_norm);
93
94 possion_norm = sqrt(possion_norm);
95 elec_continuty_norm = sqrt(elec_continuty_norm);
96 hole_continuty_norm = sqrt(hole_continuty_norm);
97 }
98
99
100 /* ----------------------------------------------------------------------------
101 * Convergence Test function for line search method.
102 */
103 PetscErrorCode LineSearch_ConvergenceTest_Mix1(SNES snes,PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm,
104 SNESConvergedReason *reason,void *dummy)
105 {
106 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)dummy;
107
108 *reason = SNES_CONVERGED_ITERATING;
109 if (!it)
110 {
111 snes->ttol = fnorm*snes->rtol;
112 gss_log.string_buf()<<" "<<"its\t"<<"| Eq(V) | "<<"| Eq(n) | "<<"| Eq(p) | "<<"|delta x|\n"
113 <<"----------------------------------------------------------------------\n";
114 gss_log.record();
115 }
116
117 gss_log.string_buf().precision(3);
118 gss_log.string_buf()<<" "<<it<<"\t"
119 <<ps->possion_norm<<" "
120 <<ps->elec_continuty_norm<<" "
121 <<ps->hole_continuty_norm<<" "
122 <<pnorm<<"\n";
123 gss_log.record();
124 gss_log.string_buf().precision(6);
125
126 if (fnorm != fnorm)
127 {
128 *reason = SNES_DIVERGED_FNORM_NAN;
129 }
130 else if (ps->possion_norm < ps->possion_abs_toler &&
131 ps->elec_continuty_norm < ps->elec_continuty_abs_toler &&
132 ps->hole_continuty_norm < ps->hole_continuty_abs_toler )
133 {
134 *reason = SNES_CONVERGED_FNORM_ABS;
135 }
136 else if (snes->nfuncs >= snes->max_funcs)
137 {
138 *reason = SNES_DIVERGED_FUNCTION_COUNT;
139 }
140
141 if (it && !*reason)
142 {
143 if (fnorm <= snes->ttol)
144 {
145 *reason = SNES_CONVERGED_FNORM_RELATIVE;
146 }
147 else if (pnorm < ps->relative_toler &&
148 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
149 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
150 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler )
151 {
152 *reason = SNES_CONVERGED_PNORM_RELATIVE;
153 }
154 }
155
156 if(*reason>0)
157 {
158 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
159 <<" "<<SNESConvergedReasons[*reason]<<" *residual norm = "<<fnorm<<"\n\n\n";
160 gss_log.record();
161 }
162
163
164 return(0);
165 }
166
167
168 /* ----------------------------------------------------------------------------
169 * Convergence Test function for trust region method.
170 */
171 PetscErrorCode TrustRegion_ConvergenceTest_Mix1(SNES snes,PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm,
172 SNESConvergedReason *reason,void *dummy)
173 {
174 SNES_TR *neP = (SNES_TR *)snes->data;
175 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)dummy;
176
177 if (!ps->its)
178 {
179 gss_log.string_buf()<<" "<<"its\t"<<"| Eq(V) | "<<"| Eq(n) | "<<"| Eq(p) | "<<"|delta x|\n"
180 <<"----------------------------------------------------------------------\n";
181 gss_log.record();
182 }
183
184 gss_log.string_buf().precision(3);
185 gss_log.string_buf()<<" "<<ps->its++<<"\t"
186 <<ps->possion_norm<<" "
187 <<ps->elec_continuty_norm<<" "
188 <<ps->hole_continuty_norm<<" "
189 <<pnorm<<"\n";
190 gss_log.record();
191 gss_log.string_buf().precision(6);
192
193 *reason = SNES_CONVERGED_ITERATING;
194
195 if (!it)
196 {
197 snes->ttol = fnorm*snes->rtol;
198 }
199
200 if (fnorm != fnorm)
201 {
202 *reason = SNES_DIVERGED_FNORM_NAN;
203 }
204 else if (neP->delta < xnorm * snes->deltatol)
205 {
206 if( ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
207 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
208 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler )
209 {
210 *reason = SNES_CONVERGED_TR_DELTA;
211 }
212 else
213 {
214 *reason = SNES_DIVERGED_LOCAL_MIN;
215 }
216 }
217 else if (ps->possion_norm < ps->possion_abs_toler &&
218 ps->elec_continuty_norm < ps->elec_continuty_abs_toler &&
219 ps->hole_continuty_norm < ps->hole_continuty_abs_toler )
220 {
221 *reason = SNES_CONVERGED_FNORM_ABS;
222 }
223 else if (snes->nfuncs >= snes->max_funcs)
224 {
225 *reason = SNES_DIVERGED_FUNCTION_COUNT;
226 }
227
228 if (it && !*reason)
229 {
230 if (fnorm <= snes->ttol)
231 {
232 *reason = SNES_CONVERGED_FNORM_RELATIVE;
233 }
234 else if (pnorm < ps->relative_toler &&
235 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
236 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
237 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler )
238 {
239 *reason = SNES_CONVERGED_PNORM_RELATIVE;
240 }
241 }
242
243 if (snes->nfuncs >= snes->max_funcs)
244 {
245 *reason = SNES_DIVERGED_FUNCTION_COUNT;
246 }
247 if(*reason>0)
248 {
249 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
250 <<" "<<SNESConvergedReasons[*reason]<<" *residual norm = "<<fnorm<<"\n\n\n";
251 gss_log.record();
252 }
253
254 return(0);
255 }
256
257
258
259 /* ----------------------------------------------------------------------------
260 * form_function_pn: This function setup DDM equation F(x)=0
261 */
262 void DDM_Mix_Solver_L1E::form_function_pn_Mix1(PetscScalar *x,PetscScalar *f)
263 {
264 // compute flux along triangle edges. semiconductor zone only
265 for(int z=0;z<zone_num;z++)
266 if(zonedata[z]->material_type==Semiconductor)
267 {
268 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
269 for(int i=0;i<zone[z].datri.size();i++)
270 {
271 Tri *ptri = &zone[z].datri[i];
272 pzonedata->F1E_Tri_ddm(ptri,x,f,zofs);
273 }
274 }
275
276 for(int z=0;z<zone_num;z++)
277 {
278 //semiconductor zone
279 if(zonedata[z]->material_type==Semiconductor)
280 {
281 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
282 // compute electrode current.
283 for(int i=0;i<pzonedata->electrode.size();i++)
284 {
285 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
286 pzonedata->F1E_mix_om_electrode_current(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
287 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
288 pzonedata->F1E_mix_stk_electrode_current(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
289 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
290 pzonedata->F1E_mix_ins_electrode_current(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
291 }
292 // process cell variables and boundaries.
293 for(int i=0;i<zone[z].davcell.size();i++)
294 {
295 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
296 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
297 pzonedata->F1E_ddm_inner(i,x,f,ODE_formula,zofs);
298 //process the electrode bc for mix simulation.
299 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
300 pzonedata->F1E_mix_ddm_ombc(i,x,f,ODE_formula,zofs,bc);
301 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
302 pzonedata->F1E_mix_ddm_stkbc(i,x,f,ODE_formula,zofs,bc);
303 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
304 pzonedata->F1E_mix_ddm_insulator_gate(i,x,f,ODE_formula,zofs,bc);
305 // the routine for interface BC do not change, share the DDML1E routine.
306 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
307 {
308 InsulatorInterfaceBC *pbc;
309 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
310 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
311 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
312 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
313 pzonedata->F1E_ddm_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
314 }
315 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HomoInterface)
316 {
317 HomoInterfaceBC *pbc;
318 pbc = dynamic_cast<HomoInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
319 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
320 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
321 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
322 pzonedata->F1E_ddm_homojunction(i,x,f,ODE_formula,zofs,bc,pz,n_node);
323 }
324 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HeteroInterface)
325 {
326 HeteroInterfaceBC *pbc;
327 pbc = dynamic_cast<HeteroInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
328 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
329 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
330 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
331 pzonedata->F1E_ddm_heterojunction(i,x,f,ODE_formula,zofs,bc,pz,n_node);
332 }
333 else //prevent some mistakes caused by inexact bc
334 pzonedata->F1E_ddm_inner(i,x,f,ODE_formula,zofs);
335 }
336 }
337
338 // Insulator zone
339 if(zonedata[z]->material_type==Insulator)
340 {
341 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
342 for(int i=0;i<zone[z].davcell.size();i++)
343 {
344 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
345 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
346 pzonedata->F1_ddm_inner(i,x,f,ODE_formula,zofs);
347 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
348 pzonedata->F1_mix_ddm_gatebc(i,x,f,ODE_formula,zofs,bc);
349 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
350 pzonedata->F1_ddm_chargebc(i,x,f,ODE_formula,zofs,bc);
351 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
352 {
353 ElectrodeInsulatorBC *pbc;
354 pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
355 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
356 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
357 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
358 pzonedata->F1_ddm_electrode_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
359 }
360 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
361 {
362 InsulatorInterfaceBC *pbc;
363 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
364 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
365 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
366 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
367 pzonedata->F1_ddm_semiconductor_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
368 }
369 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
370 {
371 InsulatorInsulatorBC *pbc;
372 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
373 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
374 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
375 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
376 pzonedata->F1_ddm_insulator_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
377 }
378 else //prevent some mistakes caused by inexact bc
379 pzonedata->F1_ddm_inner(i,x,f,ODE_formula,zofs);
380 }
381 for(int i=0;i<pzonedata->electrode.size();i++)
382 {
383 //compute gate electrode current
384 if(bc[pzonedata->electrode[i]].BCType==GateContact)
385 pzonedata->F1_mix_gate_electrode_current(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
386 //addtional equ for charge electrode
387 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
388 pzonedata->F1_charge_electrode(i,x,f,ODE_formula,zofs,bc);
389 }
390 }
391
392 // Electrode zone
393 if(zonedata[z]->material_type==Conductor)
394 {
395 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
396 for(int i=0;i<zone[z].davcell.size();i++)
397 {
398 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
399 if(!pcell->bc_index || bc[pcell->bc_index-1].BCType==NeumannBoundary)
400 pzonedata->F1_ddm_inner(i,x,f,ODE_formula,zofs);
401 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
402 pzonedata->F1_ddm_inner(i,x,f,ODE_formula,zofs);
403 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==OhmicContact)
404 {
405 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
406 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
407 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
408 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
409 pzonedata->F1_ddm_om_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
410 }
411 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==SchottkyContact)
412 {
413 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
414 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
415 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
416 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
417 pzonedata->F1_ddm_om_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
418 }
419 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==GateContact)
420 {
421 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
422 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
423 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
424 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
425 pzonedata->F1_ddm_gate_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
426 }
427 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
428 {
429 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
430 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
431 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
432 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
433 pzonedata->F1_ddm_charge_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
434 }
435 else //prevent some mistakes caused by inexact bc
436 pzonedata->F1_ddm_inner(i,x,f,ODE_formula,zofs);
437 }
438 }
439
440 }
441 }
442
443
444 /* ----------------------------------------------------------------------------
445 * form_jacobian_pn: This function setup jacobian matrix F'(x) of DDM equation F(x)=0
446 * dual carrier edition
447 */
448 void DDM_Mix_Solver_L1E::form_jacobian_pn_Mix1(PetscScalar *x, Mat *jac, Mat *jtmp)
449 {
450
451 //Compute Jacobian entries for flux along triangle edges.
452 for(int z=0;z<zone_num;z++)
453 {
454 if(zonedata[z]->material_type==Semiconductor)
455 {
456 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
457
458 // compute flux Jacobian along triangle edges
459 for(int i=0;i<zone[z].datri.size();i++)
460 {
461 Tri *ptri = &zone[z].datri[i];
462 pzonedata->J1E_Tri_ddm(ptri,x,jtmp,zofs);
463 }
464 }
465 }
466
467 MatAssemblyBegin(*jtmp,MAT_FINAL_ASSEMBLY);
468 MatAssemblyEnd(*jtmp,MAT_FINAL_ASSEMBLY);
469
470 //Compute Jacobian entries and insert into matrix.
471 for(int z=0;z<zone_num;z++)
472 {
473 if(zonedata[z]->material_type==Semiconductor)
474 {
475 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
476
477 // process cell variables and boundaries.
478 for(int i=0;i<zone[z].davcell.size();i++)
479 {
480 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
481 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
482 pzonedata->J1E_ddm_inner(i,x,jac,&JTmp,ODE_formula,zofs);
483 //process the electrode bc for mix simulation.
484 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
485 pzonedata->J1E_mix_ddm_ombc(i,x,jac,&JTmp,ODE_formula,zofs,bc);
486 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
487 pzonedata->J1E_mix_ddm_stkbc(i,x,jac,&JTmp,ODE_formula,zofs,bc);
488 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
489 pzonedata->J1E_mix_ddm_insulator_gate(i,x,jac,&JTmp,ODE_formula,zofs,bc);
490 // the routine for interface BC do not change, share the DDML1E routine.
491 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
492 {
493 InsulatorInterfaceBC *pbc;
494 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
495 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
496 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
497 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
498 pzonedata->J1E_ddm_interface(i,x,jac,&JTmp,ODE_formula,zofs,bc,pz,n_node);
499 }
500 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HomoInterface)
501 {
502 HomoInterfaceBC *pbc;
503 pbc = dynamic_cast<HomoInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
504 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
505 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
506 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
507 pzonedata->J1E_ddm_homojunction(i,x,jac,&JTmp,ODE_formula,zofs,bc,pz,n_node);
508 }
509 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==HeteroInterface)
510 {
511 HeteroInterfaceBC *pbc;
512 pbc = dynamic_cast<HeteroInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
513 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
514 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
515 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
516 pzonedata->J1E_ddm_heterojunction(i,x,jac,&JTmp,ODE_formula,zofs,bc,pz,n_node);
517 }
518 else //prevent some mistakes caused by inexact bc
519 pzonedata->J1E_ddm_inner(i,x,jac,&JTmp,ODE_formula,zofs);
520 }
521 }
522
523 // Insulator zone
524 if(zonedata[z]->material_type==Insulator)
525 {
526 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
527
528 for(int i=0;i<zone[z].davcell.size();i++)
529 {
530
531 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
532 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
533 pzonedata->J1_ddm_inner(i,x,jac,ODE_formula,zofs);
534 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
535 pzonedata->J1_mix_ddm_gatebc(i,x,jac,ODE_formula,zofs,bc);
536 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==ChargedContact)
537 pzonedata->J1_ddm_chargebc(i,x,jac,ODE_formula,zofs,bc);
538 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
539 {
540 ElectrodeInsulatorBC *pbc;
541 pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
542 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
543 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
544 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
545 pzonedata->J1_ddm_electrode_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
546 }
547 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
548 {
549 InsulatorInterfaceBC *pbc;
550 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
551 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
552 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
553 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
554 pzonedata->J1_ddm_semiconductor_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
555 }
556 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
557 {
558 InsulatorInsulatorBC *pbc;
559 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
560 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
561 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
562 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
563 pzonedata->J1_ddm_insulator_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
564 }
565 else //prevent some mistakes caused by inexact bc
566 pzonedata->J1_ddm_inner(i,x,jac,ODE_formula,zofs);
567 }
568 for(int i=0;i<pzonedata->electrode.size();i++)
569 {
570 //addtional equ for charge electrode
571 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
572 pzonedata->J1_charge_electrode(i,x,jac,ODE_formula,zofs,bc);
573 }
574 }
575
576 // Electrode zone
577 if(zonedata[z]->material_type==Conductor)
578 {
579 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
580 for(int i=0;i<zone[z].davcell.size();i++)
581 {
582 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
583 if(!pcell->bc_index || bc[pcell->bc_index-1].BCType==NeumannBoundary)
584 pzonedata->J1_ddm_inner(i,x,jac,ODE_formula,zofs);
585 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
586 pzonedata->J1_ddm_inner(i,x,jac,ODE_formula,zofs);
587 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==OhmicContact)
588 {
589 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
590 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
591 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
592 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
593 pzonedata->J1_ddm_om_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
594 }
595 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==SchottkyContact)
596 {
597 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
598 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
599 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
600 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
601 pzonedata->J1_ddm_stk_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
602 }
603 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==GateContact)
604 {
605 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
606 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
607 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
608 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
609 pzonedata->J1_ddm_gate_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
610 }
611 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
612 {
613 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
614 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
615 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
616 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
617 pzonedata->J1_ddm_charge_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
618 }
619 else //prevent some mistakes caused by inexact bc
620 pzonedata->J1_ddm_inner(i,x,jac,ODE_formula,zofs);
621 }
622 }
623
624 }
625 }
626
627
628 /* ----------------------------------------------------------------------------
629 * SNES_form_function_pn_Mix1: wrap function for petsc nonlinear solver
630 */
631 PetscErrorCode SNES_form_function_pn_Mix1(SNES snes, Vec x,Vec f,void *dummy)
632 {
633 PetscScalar *xx,*ff;
634 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)dummy;
635 VecZeroEntries(f);
636 //Get pointers to vector data.
637 VecGetArray(x,&xx);
638 VecGetArray(f,&ff);
639
640 ps->form_function_pn_Mix1(xx,ff);
641 ps->error_norm_pn_Mix1(xx,ff);
642
643 //Restore vectors
644 VecRestoreArray(x,&xx);
645 VecRestoreArray(f,&ff);
646 VecNorm(f,NORM_2,&ps->norm);
647 return 0;
648 }
649
650
651 /* ----------------------------------------------------------------------------
652 * SNES_form_jacobian_pn_Mix1: wrap function for petsc nonlinear solver
653 */
654 PetscErrorCode SNES_form_jacobian_pn_Mix1(SNES snes, Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
655 {
656 PetscScalar *xx;
657 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)dummy;
658 //Get pointer to vector data
659 VecGetArray(x,&xx);
660 //clear old matrix
661 MatZeroEntries(*jac);
662 MatZeroEntries(ps->JTmp);
663 //build matrix here
664 ps->form_jacobian_pn_Mix1(xx,jac,&ps->JTmp);
665 *flag = SAME_NONZERO_PATTERN;
666 //Restore vector
667 VecRestoreArray(x,&xx);
668 //Assemble matrix
669 MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
670 MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
671 //MatView(*jac,PETSC_VIEWER_DRAW_WORLD);
672 return 0;
673 }
674
675 PetscErrorCode LimitorByBankRose_Mix1(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
676 {
677 int it;
678 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)checkctx;
679 PetscScalar *xx;
680 PetscScalar *yy;
681 PetscScalar *ww;
682 VecGetArray(x,&xx);
683 VecGetArray(y,&yy);
684 VecGetArray(w,&ww);
685
686 SNESGetIterationNumber(snes,&it);
687 static PetscScalar K;
688 if(it==0) K=0;
689
690 Vec gk0,gk1;
691 int flag=0;
692 VecDuplicate(x,&gk0);
693 VecDuplicate(x,&gk1);
694 SNESComputeFunction(snes,x,gk0);
695 PetscScalar gk0_norm;
696 VecNorm(gk0,NORM_2,&gk0_norm);
697 while(1)
698 {
699 PetscScalar tk=1.0/(1+K*gk0_norm);
700 if(tk<1e-4) break;
701 //w=x-tk*y;
702 flag=1;
703 VecWAXPY(w,-tk,y,x);
704 SNESComputeFunction(snes,w,gk1);
705 PetscScalar gk1_norm;
706 VecNorm(gk1,NORM_2,&gk1_norm);
707 if((1-gk1_norm/gk0_norm)<0.9*tk)
708 {
709 if(K==0) K=1;
710 else K*=10;
711 }
712 else
713 {
714 K/=10;
715 break;
716 }
717 }
718
719 //prevent negative carrier density
720 int offset=0;
721 for(int z=0;z<ps->zone_num;z++)
722 {
723 if(ps->zonedata[z]->material_type==Semiconductor)
724 {
725 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
726 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
727 {
728 if(fabs(yy[offset])>1.0) {ww[offset] = xx[offset] - dsign(yy[offset])*1.0;flag=1;}
729 if (ww[offset+1] <0 ) {ww[offset+1]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
730 if (ww[offset+2] <0 ) {ww[offset+2]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
731 offset += 3;
732 }
733 }
734 else if(ps->zonedata[z]->material_type==Insulator)
735 {
736 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
737 offset += pzonedata->pzone->davcell.size();
738 offset += pzonedata->electrode.size();
739 }
740 else if(ps->zonedata[z]->material_type==Conductor)
741 {
742 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
743 offset += pzonedata->pzone->davcell.size();
744 }
745 }
746
747 if(flag)
748 {
749 *changed_y = PETSC_FALSE;
750 *changed_w = PETSC_TRUE;
751 }
752 VecDestroy(gk0);
753 VecDestroy(gk1);
754 VecRestoreArray(x,&xx);
755 VecRestoreArray(y,&yy);
756 VecRestoreArray(w,&ww);
757 return(0);
758 }
759
760
761 PetscErrorCode LimitorByPotential_Mix1(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
762 {
763 PetscScalar *xx;
764 PetscScalar *yy;
765 PetscScalar *ww;
766 VecGetArray(x,&xx);
767 VecGetArray(y,&yy);
768 VecGetArray(w,&ww);
769 PetscScalar dV_max=0;
770
771 //search for dV_max;
772 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)checkctx;
773 int offset=0;
774 for(int z=0;z<ps->zone_num;z++)
775 {
776 if(ps->zonedata[z]->material_type==Semiconductor)
777 {
778 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
779 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
780 {
781 if(fabs(yy[offset])>dV_max) {dV_max=fabs(yy[offset]); }
782 offset += 3;
783 }
784 }
785 else if(ps->zonedata[z]->material_type==Insulator)
786 {
787 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
788 offset += pzonedata->pzone->davcell.size();
789 offset += pzonedata->electrode.size();
790 }
791 else if(ps->zonedata[z]->material_type==Conductor)
792 {
793 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
794 offset += pzonedata->pzone->davcell.size();
795 }
796 }
797 if(dV_max<1e-2) return(0);
798
799 //do logarithmic potential damping;
800 PetscScalar Vt=0.026*ps->LatticeTemp;
801 PetscScalar f = log(1+dV_max/Vt)/(dV_max/Vt);
802
803 //prevent negative carrier density
804 offset=0;
805 for(int z=0;z<ps->zone_num;z++)
806 {
807 if(ps->zonedata[z]->material_type==Semiconductor)
808 {
809 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
810 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
811 {
812 ww[offset] = xx[offset] - f*yy[offset];
813 if (ww[offset+1] <0 ) ww[offset+1]=1.0*pow(ps->scale_unit.s_centimeter,-3);
814 if (ww[offset+2] <0 ) ww[offset+2]=1.0*pow(ps->scale_unit.s_centimeter,-3);
815 offset += 3;
816 }
817 }
818 else if(ps->zonedata[z]->material_type==Insulator)
819 {
820 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
821 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
822 {
823 ww[offset] = xx[offset] - f*yy[offset];
824 offset += 1;
825 }
826 offset += pzonedata->electrode.size();
827 }
828 else if(ps->zonedata[z]->material_type==Conductor)
829 {
830 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
831 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
832 {
833 ww[offset] = xx[offset] - f*yy[offset];
834 offset += 1;
835 }
836 }
837 }
838
839 VecRestoreArray(x,&xx);
840 VecRestoreArray(y,&yy);
841 VecRestoreArray(w,&ww);
842
843 *changed_y = PETSC_FALSE;
844 *changed_w = PETSC_TRUE;
845 return(0);
846 }
847
848
849 PetscErrorCode LimitorNonNegativeCarrier_Mix1(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
850 {
851 PetscScalar *xx;
852 PetscScalar *yy;
853 PetscScalar *ww;
854 VecGetArray(x,&xx);
855 VecGetArray(y,&yy);
856 VecGetArray(w,&ww);
857
858 DDM_Mix_Solver_L1E *ps = (DDM_Mix_Solver_L1E*)checkctx;
859 int flag=0;
860 int offset=0;
861 for(int z=0;z<ps->zone_num;z++)
862 {
863 if(ps->zonedata[z]->material_type==Semiconductor)
864 {
865 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
866 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
867 {
868 if(fabs(yy[offset])>1.0) {ww[offset] = xx[offset] - dsign(yy[offset])*1.0;flag=1;}
869 if (ww[offset+1] <0 ) {ww[offset+1]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
870 if (ww[offset+2] <0 ) {ww[offset+2]=1.0*pow(ps->scale_unit.s_centimeter,-3);flag=1;}
871 offset += 3;
872 }
873 }
874 else if(ps->zonedata[z]->material_type==Insulator)
875 {
876 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
877 offset += pzonedata->pzone->davcell.size();
878 offset += pzonedata->electrode.size();
879 }
880 else if(ps->zonedata[z]->material_type==Conductor)
881 {
882 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
883 offset += pzonedata->pzone->davcell.size();
884 }
885 }
886 VecRestoreArray(x,&xx);
887 VecRestoreArray(y,&yy);
888 VecRestoreArray(w,&ww);
889
890 if(flag)
891 {
892 *changed_y = PETSC_FALSE;
893 *changed_w = PETSC_TRUE;
894 }
895
896 return(0);
897 }
898
899
900 /* ----------------------------------------------------------------------------
901 * DDM_Mix_Solver_L1E::init_solver: This function do initial setup for nonlinear solver
902 */
903 int DDM_Mix_Solver_L1E::init_solver(SolveDefine &sv)
904 {
905 gss_log.string_buf()<<"DDM Level 1E Device/Circuit Mix Solver init...\n";
906 gss_log.record();
907 psv = &sv;
908 //connect to ngspice
909 if(NDEVServer(sv.port,listener,client))
910 {
911 gss_log.string_buf()<<"I can't connect to ngspice, mix simulation stopped.\n";
912 gss_log.record();
913 return 1;
914 }
915 //get ngspice terminal information
916 recv(client,&(Deviceinfo),sizeof(sDeviceinfo),MSG_FLAG);
917 gss_log.string_buf()<<"Ngspice demands device "<<Deviceinfo.NDEVname<<" to be calculated.\n";
918 gss_log.record();
919 for(int i=0;i<Deviceinfo.term;i++)
920 {
921 recv(client,&(PINinfos[i]),sizeof(sPINinfo),MSG_FLAG);
922 if(!bc.Get_electrode_pointer_nocase(PINinfos[i].name))
923 {
924 gss_log.string_buf()<<"I can't link terminal name "<<PINinfos[i].name<<" to GSS eletrode boundary.\n";
925 gss_log.record();
926 return 1;
927 }
928 gss_log.string_buf()<<" Terminal "<<i<<": "<<PINinfos[i].name<<"\n";
929 gss_log.record();
930 }
931
932 //set Tolerances
933 relative_toler = sv.relative_toler;
934 toler_relax = sv.toler_relax;
935 possion_abs_toler = sv.possion_abs_toler;
936 elec_continuty_abs_toler = sv.elec_continuty_abs_toler;
937 hole_continuty_abs_toler = sv.hole_continuty_abs_toler;
938
939 // compute the scale of problem.
940 zofs.resize(zone_num+1);
941 zofs[0] = 0;
942 for(int i=0;i<zone_num;i++)
943 {
944 if(zonedata[i]->material_type==Semiconductor) //semiconductor zone
945 {
946 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
947 N += 3*pzonedata->pzone->davcell.size();
948 //NO extra electrode equation needs.
949 zofs[i+1] = N;
950 }
951 else if(zonedata[i]->material_type==Insulator) //Insulator zone
952 {
953 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
954 N += zone[i].davcell.size();
955 N += pzonedata->electrode.size(); //additional equation for charged electrode
956 zofs[i+1] = N;
957 }
958 else if(zonedata[i]->material_type==Conductor) //Electrode zone
959 {
960 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
961 N += zone[i].davcell.size();
962 zofs[i+1] = N;
963 }
964 else //other zones, such as PML and vacuum
965 {
966 zofs[i+1] = N;
967 }
968 }
969
970 VecCreateSeq(PETSC_COMM_SELF,N,&x);
971 VecDuplicate(x,&r);
972 for(int i=0;i<Deviceinfo.term;i++)
973 {
974 VecDuplicate(x,&PINconds[i].pdI_pdw);
975 VecDuplicate(x,&PINconds[i].pdF_pdV);
976 VecDuplicate(x,&PINconds[i].pdw_pdV);
977 }
978
979 //Evaluate initial guess,importnat for newton solver
980 PetscScalar init_value[3];
981 int index[3];
982 int offset=0;
983 for(int z=0;z<zone_num;z++)
984 {
985 if(zonedata[z]->material_type==Semiconductor)
986 {
987 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
988 for(int i=0;i<zone[z].davcell.size();i++)
989 {
990 index[0] = offset+0;
991 index[1] = offset+1;
992 index[2] = offset+2;
993 init_value[0] = pzonedata->fs[i].P;
994 init_value[1] = pzonedata->fs[i].n;
995 init_value[2] = pzonedata->fs[i].p;
996 VecSetValues(x,3,index,init_value,INSERT_VALUES);
997 offset += 3;
998 }
999 }
1000
1001 if(zonedata[z]->material_type==Insulator)
1002 {
1003 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1004 for(int i=0;i<zone[z].davcell.size();i++)
1005 {
1006 VecSetValue(x,offset,pzonedata->fs[i].P,INSERT_VALUES);
1007 offset += 1;
1008 }
1009 for(int j=0;j<pzonedata->electrode.size();j++)
1010 {
1011 VecSetValue(x,offset,0,INSERT_VALUES);
1012 offset += 1;
1013 }
1014 }
1015
1016 if(zonedata[z]->material_type==Conductor)
1017 {
1018 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1019 for(int i=0;i<zone[z].davcell.size();i++)
1020 {
1021 VecSetValue(x,offset,pzonedata->fs[i].P,INSERT_VALUES);
1022 offset += 1;
1023 }
1024 }
1025
1026 }
1027 //
1028 VecAssemblyBegin(x);
1029 VecAssemblyEnd(x);
1030
1031 // Create Jacobian matrix data structure
1032 // pre-alloc approximate memory
1033 int *nnz = new int[N];
1034 int *p = nnz;
1035 for(int i=0;i<zone_num;i++)
1036 {
1037 //semiconductor zone
1038 if(zonedata[i]->material_type==Semiconductor)
1039 {
1040 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1041 for(int j=0;j<pzonedata->pzone->davcell.size();j++)
1042 {
1043 const VoronoiCell* pcell = pzonedata->pzone->davcell.GetPointer(j);
1044 //inner node, alloc exact memory.
1045 if(!pcell->bc_index || bc[pcell->bc_index-1].psegment->interface==-1)
1046 {
1047 *p++ = 3*(pzonedata->pzone->davcell[j].nb_num+1);
1048 *p++ = 3*(pzonedata->pzone->davcell[j].nb_num+1);
1049 *p++ = 3*(pzonedata->pzone->davcell[j].nb_num+1);
1050 }
1051 else //interface node, slightly more than needed.
1052 {
1053 *p++ = 3*(2*pzonedata->pzone->davcell[j].nb_num+4);
1054 *p++ = 3*(2*pzonedata->pzone->davcell[j].nb_num+4);
1055 *p++ = 3*(2*pzonedata->pzone->davcell[j].nb_num+4);
1056 }
1057 }
1058 }
1059 //Insulator zones, poisson equation only
1060 if(zonedata[i]->material_type==Insulator)
1061 {
1062 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1063 for(int j=0;j<zone[i].davcell.size();j++)
1064 *p++ = 1*(zone[i].davcell[j].nb_num+1);
1065 for(int j=0;j<pzonedata->electrode.size();j++)
1066 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
1067 }
1068 //Electrode zones, poisson equation only
1069 if(zonedata[i]->material_type==Conductor)
1070 {
1071 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
1072 for(int j=0;j<zone[i].davcell.size();j++)
1073 *p++ = 1*(zone[i].davcell[j].nb_num+1);
1074 }
1075 }
1076
1077 MatCreate(PETSC_COMM_SELF,&J);
1078 MatSetSizes(J,N,N,N,N);
1079 if(sv.LS=="superlu")
1080 {
1081 #ifdef PETSC_HAVE_SUPERLU
1082 MatSetType(J,MATSUPERLU);
1083 #else
1084 MatSetType(J,MATSEQAIJ);
1085 #endif
1086 }
1087 else if(sv.LS=="umfpack")
1088 {
1089 #ifdef PETSC_HAVE_UMFPACK
1090 MatSetType(J,MATUMFPACK);
1091 #else
1092 MatSetType(J,MATSEQAIJ);
1093 #endif
1094 }
1095 else
1096 {
1097 MatSetType(J,MATSEQAIJ);
1098 }
1099 MatSeqAIJSetPreallocation(J,0,nnz);
1100
1101 MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,0,nnz,&JTmp);
1102 delete [] nnz;
1103
1104 SNESCreate(PETSC_COMM_WORLD,&snes);
1105 SNESSetFunction(snes,r,SNES_form_function_pn_Mix1,this);
1106 SNESSetJacobian(snes,J,J,SNES_form_jacobian_pn_Mix1,this);
1107
1108 // set the newton method
1109 SNESSetType(snes,SNESLS); //default method
1110
1111 // the maximum number of iterations
1112 PetscInt maxit = sv.maxit;
1113
1114 SNESSetType(snes,SNESLS); //default method
1115 if(sv.NS==LineSearch || sv.NS==Basic)
1116 {
1117 if(sv.NS==LineSearch)
1118 SNESLineSearchSet(snes,SNESLineSearchCubic,PETSC_NULL);
1119 if(sv.NS==Basic)
1120 SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);
1121
1122 if(sv.Damping==DampingBankRose)
1123 SNESLineSearchSetPostCheck(snes,LimitorByBankRose_Mix1,this);
1124 else if(sv.Damping==DampingPotential)
1125 SNESLineSearchSetPostCheck(snes,LimitorByPotential_Mix1,this);
1126 else
1127 SNESLineSearchSetPostCheck(snes,LimitorNonNegativeCarrier_Mix1,this);
1128
1129 SNESSetTolerances(snes,1e-12*N,1e-14,1e-9,maxit,1000);
1130 SNESSetConvergenceTest(snes,LineSearch_ConvergenceTest_Mix1,this);
1131 }
1132 else if(sv.NS==TrustRegion)
1133 {
1134 SNESSetType(snes,SNESTR);
1135 SNESSetTolerances(snes,1e-12*N,1e-14,1e-9,maxit,1000);
1136 SNESSetTrustRegionTolerance(snes,1e-30);
1137 SNESSetConvergenceTest(snes,TrustRegion_ConvergenceTest_Mix1,this);
1138 //must set TR delta0 to a sufficient large value, or TR can't find real solution.
1139 SNES_TR *neP = (SNES_TR*)snes->data;
1140 neP->delta0 = N;
1141 }
1142 SNESGetKSP(snes,&ksp);
1143 KSPGetPC(ksp,&pc);
1144
1145 if(sv.LS=="lu"||sv.LS=="superlu"||sv.LS=="umfpack")
1146 {
1147 KSPSetType(ksp,KSPPREONLY);
1148 PCSetType(pc,PCLU);
1149 PCFactorSetReuseOrdering(pc,PETSC_TRUE);
1150 PCFactorSetReuseFill(pc,PETSC_TRUE);
1151 PCFactorSetPivoting(pc,1.0);
1152 PCFactorReorderForNonzeroDiagonal(pc,1e-14);
1153 PCFactorSetShiftNonzero(pc,1e-12);
1154 }
1155 else
1156 {
1157 KSPSetType(ksp,sv.LS.c_str());
1158 if(sv.LS=="gmres") KSPGMRESSetRestart(ksp,150);
1159 KSPSetTolerances(ksp,1e-14*N,1e-25*N,PETSC_DEFAULT,max(35,N/10));
1160 PCSetType(pc,PCILU);
1161 PCFactorSetLevels(pc,3);
1162 PCFactorSetShiftNonzero(pc,1e-14);
1163 PCFactorReorderForNonzeroDiagonal(pc,1e-12);
1164 }
1165
1166 KSPSetFromOptions(ksp);
1167 SNESSetFromOptions(snes);
1168 gss_log.string_buf()<<"DDM Level 1E Device/Circuit Mix Solver init ok.\n";
1169 gss_log.record();
1170 return 0;
1171 }
1172
1173
1174 /* ----------------------------------------------------------------------------
1175 * DDM_Mix_Solver_L1E::do_solve: This function receive ngspice device information.
1176 */
1177 int DDM_Mix_Solver_L1E::do_solve(SolveDefine &sv)
1178 {
1179 ODE_formula.clock = -1e100;
1180 int bytes;
1181 for(;;)
1182 {
1183 if (sigsetjmp(net_buf, 1) == 1) break;
1184 //receive solver information
1185 bytes=recv(client,&(CKTInfo),sizeof(sCKTinfo),MSG_FLAG);
1186 //if ngspice terminated, stop loop.
1187 if(bytes<sizeof(sCKTinfo))
1188 {
1189 gss_log.string_buf()<<"Lost connection to remote ngspice.\n";
1190 gss_log.record();
1191 break;
1192 }
1193 switch(CKTInfo.DEV_CALL)
1194 {
1195 case NDEV_LOAD : DEV_LOAD(); break;
1196 case NDEV_ACCEPT : DEV_ACCEPT(); break;
1197 case NDEV_CONVERGINCE_TEST : DEV_CONV_TEST(); break;
1198 default : break;
1199 }
1200 }
1201 return 0;
1202 }
1203
1204
1205 int DDM_Mix_Solver_L1E::DEV_LOAD()
1206 {
1207 PetscScalar *xx;
1208 //receive terminal voltage of device
1209 for(int i=0;i<Deviceinfo.term;i++)
1210 {
1211 recv(client,&PINinfos[i],sizeof(sPINinfo),MSG_FLAG);
1212 bc.Set_Vapp_nocase(PINinfos[i].name,PINinfos[i].V);
1213 }
1214
1215 //transient calculation settings
1216 if(CKTInfo.CKTmode & MODETRAN)
1217 {
1218 ODE_formula.TimeDependent = true;
1219 ODE_formula.BDF_Type = BDF2;
1220 if(CKTInfo.CKTmode & MODEINITTRAN) //the first step of transient simulation
1221 {
1222 gss_log.string_buf()<<"----------------------------------------------------------------------\n";
1223 gss_log.string_buf()<<"NGSPICE Transient Mode Start Time = "<<CKTInfo.time<<"s dt = "<<CKTInfo.dt<<"s\n";
1224 for(int i=0;i<Deviceinfo.term;i++)
1225 {
1226 gss_log.string_buf()<<"Set V("<<PINinfos[i].name<<") = "<<PINinfos[i].V<<"V\n";
1227 gss_log.record();
1228 }
1229 gss_log.string_buf()<<"----------------------------------------------------------------------\n\n";
1230 gss_log.record();
1231 ODE_formula.BDF2_restart = true;
1232 }
1233 else if(CKTInfo.CKTmode & MODEINITPRED) //the time matching flag
1234 {
1235 gss_log.string_buf()<<"----------------------------------------------------------------------\n";
1236 gss_log.string_buf()<<"NGSPICE Transient Mode Update Time = "<<CKTInfo.time<<"s dt = "<<CKTInfo.dt<<"s\n";
1237 for(int i=0;i<Deviceinfo.term;i++)
1238 {
1239 gss_log.string_buf()<<"Set V("<<PINinfos[i].name<<") = "<<PINinfos[i].V<<"V\n";
1240 gss_log.record();
1241 }
1242 if(ODE_formula.clock > CKTInfo.time*scale_unit.s_second)
1243 {
1244 time_back_recovery();
1245 ODE_formula.BDF2_restart = true;
1246 gss_log.string_buf()<<"NGSPICE back trace...\n";
1247 gss_log.record();
1248 }
1249 else
1250 {
1251 ODE_formula.BDF2_restart = false;
1252 }
1253 gss_log.string_buf()<<"----------------------------------------------------------------------\n\n";
1254 gss_log.record();
1255 }
1256 ODE_formula.clock = CKTInfo.time*scale_unit.s_second;
1257 // call the real computation routine
1258 tran_solve();
1259 }
1260 //DC calculation settings
1261 else if(CKTInfo.CKTmode & MODEDC)
1262 {
1263 gss_log.string_buf()<<"----------------------------------------------------------------------\n";
1264 gss_log.string_buf()<<"NGSPICE DC Mode\n";
1265 for(int i=0;i<Deviceinfo.term;i++)
1266 {
1267 gss_log.string_buf()<<"Set V("<<PINinfos[i].name<<") = "<<PINinfos[i].V<<"V\n";
1268 gss_log.record();
1269 }
1270 gss_log.string_buf()<<"----------------------------------------------------------------------\n\n";
1271 gss_log.record();
1272 ODE_formula.TimeDependent = false;
1273 ODE_formula.dt = 1e100;
1274 // call the real computation routine
1275 dc_solve();
1276 }
1277
1278 //get pdI/pdw, pdI/pdV and pdF/pdV for each electrode
1279 VecGetArray(x,&xx);
1280 for(int i=0;i<Deviceinfo.term;i++)
1281 for(int z=0;z<zone_num;z++)
1282 {
1283 if(zonedata[z]->material_type==Semiconductor)
1284 {
1285 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1286 for(int j=0;j<pzonedata->electrode.size();j++)
1287 if(bc.is_electrode_label_nocase(pzonedata->electrode[j], PINinfos[i].name))
1288 {
1289 if(bc[pzonedata->electrode[j]].BCType==OhmicContact)
1290 pzonedata->F1E_mix_om_electrode_Load(pzonedata->electrode[j],PINinfos[i].I,PINconds[i].pdI_pdV,
1291 PINconds[i].pdI_pdw,PINconds[i].pdF_pdV,
1292 xx,&JTmp,ODE_formula,zofs,bc,DeviceDepth);
1293 if(bc[pzonedata->electrode[j]].BCType==SchottkyContact)
1294 pzonedata->F1E_mix_stk_electrode_Load(pzonedata->electrode[j],PINinfos[i].I,PINconds[i].pdI_pdV,
1295 PINconds[i].pdI_pdw,PINconds[i].pdF_pdV,
1296 xx,&JTmp,ODE_formula,zofs,bc,DeviceDepth);
1297 if(bc[pzonedata->electrode[j]].BCType==InsulatorContact)
1298 pzonedata->F1E_mix_ins_electrode_Load(pzonedata->electrode[j],PINinfos[i].I,PINconds[i].pdI_pdV,
1299 PINconds[i].pdI_pdw,PINconds[i].pdF_pdV,
1300 xx,&JTmp,ODE_formula,zofs,bc,DeviceDepth);
1301 }
1302 }
1303 //
1304 if(zonedata[z]->material_type==Insulator)
1305 {
1306 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
1307 for(int j=0;j<pzonedata->electrode.size();j++)
1308 if(bc.is_electrode_label_nocase(pzonedata->electrode[j], PINinfos[i].name))
1309 {
1310 pzonedata->F1_mix_gate_electrode_Load(pzonedata->electrode[j],PINinfos[i].I,PINconds[i].pdI_pdV,
1311 PINconds[i].pdI_pdw,PINconds[i].pdF_pdV,
1312 xx,&J,ODE_formula,zofs,bc,DeviceDepth);
1313
1314 }
1315 }
1316 }
1317 VecRestoreArray(x,&xx);
1318
1319 //save conductance and rhs current source.
1320 for(int i=0;i<Deviceinfo.term;i++)
1321 {
1322 PetscScalar Ig=0;
1323 for(int j=0;j<Deviceinfo.term;j++)
1324 {
1325 KSPSolve(ksp,PINconds[j].pdF_pdV,PINconds[i].pdw_pdV);
1326 PetscScalar pdI_pdV;
1327 VecDot(PINconds[i].pdI_pdw,PINconds[i].pdw_pdV,&pdI_pdV);
1328 PINinfos[i].dI_dV[j] = pdI_pdV;
1329 PINinfos[i].dI_dV[j] += PINconds[i].pdI_pdV;
1330 PINinfos[i].dI_dV[j] /= scale_unit.s_A; //scale the element of dI/dV.
1331 Ig += PINinfos[i].dI_dV[j]*PINinfos[j].V;
1332 }
1333 PINinfos[i].I = Ig - PINinfos[i].I/scale_unit.s_A;
1334 }
1335
1336 //sent conductance matrix and rhs current back to ngspice.
1337 for(int i=0;i<Deviceinfo.term;i++)
1338 send(client,&PINinfos[i],sizeof(sPINinfo),0);
1339
1340 return 0;
1341 }
1342
1343
1344 int DDM_Mix_Solver_L1E::DEV_ACCEPT()
1345 {
1346 solution_update();
1347 return 0;
1348 }
1349
1350
1351 int DDM_Mix_Solver_L1E::DEV_CONV_TEST()
1352 {
1353 CKTInfo.convergence_flag = reason;
1354 send(client,&(CKTInfo),sizeof(sCKTinfo),0);
1355 return 0;
1356 }
1357
1358
1359 int DDM_Mix_Solver_L1E::tran_solve()
1360 {
1361 PetscInt rework = 1;
1362 PetscInt max_rework = 64;
1363 PetscInt Converged = 0;
1364
1365 for(int i=0;i<zone_num;i++)
1366 if(zonedata[i]->material_type == Semiconductor)
1367 {
1368 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1369 pzonedata->HighFieldMobility = psv->HighFieldMobility;
1370 pzonedata->ImpactIonization = psv->ImpactIonization;
1371 pzonedata->IIType = psv->IIType;
1372 pzonedata->BandBandTunneling = psv->BandBandTunneling;
1373 pzonedata->IncompleteIonization = psv->IncompleteIonization;
1374 pzonedata->QuantumMechanical = psv->QuantumMechanical;
1375 pzonedata->Fermi = psv->Fermi;
1376 pzonedata->EJModel = psv->EJModel;
1377 }
1378
1379 do
1380 {
1381 reason = SNES_CONVERGED_ITERATING;
1382 for(int w=1;w<=rework;w++)
1383 {
1384 ODE_formula.dt = w*CKTInfo.dt*scale_unit.s_second/rework;
1385
1386 for(int i=0;i<Deviceinfo.term;i++)
1387 {
1388 PetscScalar V_step = PINinfos[i].V_old+(PINinfos[i].V-PINinfos[i].V_old)*w/rework;
1389 bc.Set_Vapp_nocase(PINinfos[i].name,V_step);
1390 }
1391
1392 SNESSolve(snes,PETSC_NULL,x);
1393 SNESGetConvergedReason(snes,&reason);
1394 if(reason<0)
1395 {
1396 gss_log.string_buf()<<"------> GSS mixed solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1397 gss_log.record();
1398 diverged_recovery();
1399 rework*=2;
1400 break;
1401 }
1402 }
1403 if(reason>0) {Converged = 1; break;}
1404 if(rework>max_rework) {Converged = 0; break;}
1405 } while(1);
1406
1407 return Converged;
1408 }
1409
1410
1411 int DDM_Mix_Solver_L1E::dc_solve()
1412 {
1413 PetscInt rework = 1;
1414 PetscInt max_rework = 64;
1415 PetscInt Converged = 0;
1416
1417 for(int i=0;i<zone_num;i++)
1418 if(zonedata[i]->material_type == Semiconductor)
1419 {
1420 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1421 pzonedata->HighFieldMobility = psv->HighFieldMobility;
1422 pzonedata->ImpactIonization = psv->ImpactIonization;
1423 pzonedata->IIType = psv->IIType;
1424 pzonedata->BandBandTunneling = psv->BandBandTunneling;
1425 pzonedata->IncompleteIonization = psv->IncompleteIonization;
1426 pzonedata->QuantumMechanical = psv->QuantumMechanical;
1427 pzonedata->Fermi = psv->Fermi;
1428 pzonedata->EJModel = psv->EJModel;
1429 }
1430
1431 do
1432 {
1433 reason = SNES_CONVERGED_ITERATING;
1434 for(int w=1;w<=rework;w++)
1435 {
1436 for(int i=0;i<Deviceinfo.term;i++)
1437 {
1438 PetscScalar V_step = PINinfos[i].V_old+(PINinfos[i].V-PINinfos[i].V_old)*w/rework;
1439 bc.Set_Vapp_nocase(PINinfos[i].name,V_step);
1440 }
1441
1442 SNESSolve(snes,PETSC_NULL,x);
1443 SNESGetConvergedReason(snes,&reason);
1444 if(reason<0)
1445 {
1446 gss_log.string_buf()<<"------> GSS mixed solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1447 gss_log.record();
1448 diverged_recovery();
1449 rework*=2;
1450 break;
1451 }
1452 }
1453 if(reason>0) {Converged = 1; break;}
1454 if(rework>max_rework) {Converged = 0; break;}
1455 } while(1);
1456
1457 return Converged;
1458 }
1459
1460
1461 /* ----------------------------------------------------------------------------
1462 * DDM_Mix_Solver_L1E::solution_update: This function restore solution data from SNES
1463 */
1464 void DDM_Mix_Solver_L1E::solution_update()
1465 {
1466 //------------------------------------------------------------
1467 PetscScalar *xx;
1468 VecGetArray(x,&xx);
1469 int offset=0;
1470 for(int z=0;z<zone_num;z++)
1471 {
1472 if(zonedata[z]->material_type==Semiconductor)
1473 {
1474 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1475 PetscScalar kb = pzonedata->mt->kb;
1476 PetscScalar e = pzonedata->mt->e;
1477 for(int i=0;i<zone[z].davcell.size();i++)
1478 {
1479 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
1480 pzonedata->fs[i].n_last = pzonedata->fs[i].n;
1481 pzonedata->fs[i].p_last = pzonedata->fs[i].p;
1482
1483 pzonedata->fs[i].P = xx[offset+0];
1484 pzonedata->fs[i].n = xx[offset+1];
1485 pzonedata->fs[i].p = xx[offset+2];
1486
1487 pzonedata->mt->mapping(&pzonedata->pzone->danode[i],&pzonedata->aux[i],0);
1488 PetscScalar nie = pzonedata->mt->band->nie(pzonedata->fs[i].T);
1489 pzonedata->aux[i].phi_intrinsic = pzonedata->fs[i].P + pzonedata->aux[i].affinity +
1490 kb*pzonedata->fs[i].T/e*log(pzonedata->aux[i].Nc/nie);
1491 pzonedata->aux[i].phin = pzonedata->aux[i].phi_intrinsic - log(fabs(pzonedata->fs[i].n)/nie)*kb*pzonedata->fs[i].T/e;
1492 pzonedata->aux[i].phip = pzonedata->aux[i].phi_intrinsic + log(fabs(pzonedata->fs[i].p)/nie)*kb*pzonedata->fs[i].T/e;
1493
1494 offset += 3;
1495 }
1496 pzonedata->F1E_efield_update(xx,zofs,bc,zonedata);
1497 for(int i=0;i<pzonedata->electrode.size();i++)
1498 {
1499 PetscScalar I = bc.Get_pointer(pzonedata->electrode[i])->Get_Current_new();
1500 bc.Get_pointer(pzonedata->electrode[i])->Set_Current(I);
1501 }
1502 }
1503 if(zonedata[z]->material_type==Insulator)
1504 {
1505 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1506 for(int i=0;i<zone[z].davcell.size();i++)
1507 {
1508 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
1509 pzonedata->fs[i].P = xx[offset];
1510 offset += 1;
1511 }
1512 pzonedata->F1_efield_update(xx,zofs,bc,zonedata);
1513 for(int i=0;i<pzonedata->electrode.size();i++)
1514 {
1515 PetscScalar I = bc.Get_pointer(pzonedata->electrode[i])->Get_Current_new();
1516 bc.Get_pointer(pzonedata->electrode[i])->Set_Current(I);
1517 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
1518 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential(xx[offset]);
1519 offset += 1;
1520 }
1521 }
1522 if(zonedata[z]->material_type==Conductor)
1523 {
1524 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1525 for(int i=0;i<zone[z].davcell.size();i++)
1526 {
1527 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
1528 pzonedata->fs[i].P = xx[offset];
1529 offset += 1;
1530 }
1531 }
1532 }
1533 ODE_formula.dt_last = ODE_formula.dt;
1534
1535 VecRestoreArray(x,&xx);
1536 }
1537
1538
1539 /* ----------------------------------------------------------------------------
1540 * DDM_Mix_Solver_L1E::time_back_recovery: This function recovery latest solution data
1541 * if NGSPICE time drops back.
1542 */
1543 void DDM_Mix_Solver_L1E::time_back_recovery()
1544 {
1545 //------------------------------------------------------------
1546 PetscScalar *xx;
1547 VecGetArray(x,&xx);
1548 int offset=0;
1549 for(int z=0;z<zone_num;z++)
1550 {
1551 if(zonedata[z]->material_type==Semiconductor)
1552 {
1553 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1554 for(int i=0;i<zone[z].davcell.size();i++)
1555 {
1556 xx[offset+0] = pzonedata->fs[i].P;
1557 xx[offset+1] = pzonedata->fs[i].n;
1558 xx[offset+2] = pzonedata->fs[i].p;
1559 pzonedata->fs[i].P = pzonedata->fs[i].P_last;
1560 pzonedata->fs[i].n = pzonedata->fs[i].n_last;
1561 pzonedata->fs[i].p = pzonedata->fs[i].p_last;
1562 offset += 3;
1563 }
1564 }
1565 if(zonedata[z]->material_type==Insulator)
1566 {
1567 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1568 for(int i=0;i<zone[z].davcell.size();i++)
1569 {
1570 xx[offset] = pzonedata->fs[i].P;
1571 pzonedata->fs[i].P = pzonedata->fs[i].P_last;
1572 offset += 1;
1573 }
1574 for(int i=0;i<pzonedata->electrode.size();i++)
1575 {
1576 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
1577 xx[offset] = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
1578 else
1579 xx[offset] = 0;
1580 offset += 1;
1581 }
1582 }
1583 if(zonedata[z]->material_type==Conductor)
1584 {
1585 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1586 for(int i=0;i<zone[z].davcell.size();i++)
1587 {
1588 xx[offset] = pzonedata->fs[i].P;
1589 pzonedata->fs[i].P = pzonedata->fs[i].P_last;
1590 offset += 1;
1591 }
1592 }
1593 }
1594 VecRestoreArray(x,&xx);
1595 }
1596
1597
1598
1599 /* ----------------------------------------------------------------------------
1600 * DDM_Mix_Solver_L1E::diverged_recovery: This function recovery latest solution data
1601 * if SNES diverged.
1602 */
1603 void DDM_Mix_Solver_L1E::diverged_recovery()
1604 {
1605 //------------------------------------------------------------
1606 PetscScalar *xx;
1607 VecGetArray(x,&xx);
1608 int offset=0;
1609 for(int z=0;z<zone_num;z++)
1610 {
1611 if(zonedata[z]->material_type==Semiconductor)
1612 {
1613 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1614 for(int i=0;i<zone[z].davcell.size();i++)
1615 {
1616 xx[offset+0] = pzonedata->fs[i].P;
1617 xx[offset+1] = pzonedata->fs[i].n;
1618 xx[offset+2] = pzonedata->fs[i].p;
1619 offset += 3;
1620 }
1621 }
1622 if(zonedata[z]->material_type==Insulator)
1623 {
1624 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1625 for(int i=0;i<zone[z].davcell.size();i++)
1626 {
1627 xx[offset] = pzonedata->fs[i].P;
1628 offset += 1;
1629 }
1630 for(int i=0;i<pzonedata->electrode.size();i++)
1631 {
1632 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
1633 xx[offset] = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
1634 else
1635 xx[offset] = 0;
1636 offset += 1;
1637 }
1638 }
1639 if(zonedata[z]->material_type==Conductor)
1640 {
1641 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1642 for(int i=0;i<zone[z].davcell.size();i++)
1643 {
1644 xx[offset] = pzonedata->fs[i].P;
1645 offset += 1;
1646 }
1647 }
1648 }
1649 VecRestoreArray(x,&xx);
1650 }
1651
1652
1653 /* ----------------------------------------------------------------------------
1654 * DDM_Mix_Solver_L1E::destroy_solver: This function do destroy the nonlinear solver
1655 */
1656 int DDM_Mix_Solver_L1E::destroy_solver(SolveDefine &sv)
1657 {
1658 // free work space
1659 N = 0;
1660 zofs.clear();
1661 for(int i=0;i<Deviceinfo.term;i++)
1662 {
1663 VecDestroy(PINconds[i].pdI_pdw);
1664 VecDestroy(PINconds[i].pdF_pdV);
1665 VecDestroy(PINconds[i].pdw_pdV);
1666 }
1667 VecDestroy(x);
1668 VecDestroy(r);
1669 MatDestroy(J);
1670 MatDestroy(JTmp);
1671 SNESDestroy(snes);
1672
1673 close(client);
1674 close(listener);
1675 return 0;
1676 }
1677
General-purpose Semiconductor Simulator