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