?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
initial checkin, based on GSS 0.46 CVS
[gss-tcad.git] / src / solver / qddm1e / qddm_nt1e.cc
bloba119429ff9e38e0cc27f83e90078049818c4b922
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: July 20, 2007 */
14 /* */
15 /* Gong Ding */
16 /* gdiso@ustc.edu */
17 /* NINT, No.69 P.O.Box, Xi'an City, China */
18 /* */
19 /*****************************************************************************/
20
21 #include "mathfunc.h"
22 #include "qddm_nt1e.h"
23 #include "private/kspimpl.h"
24 #include "private/snesimpl.h"
25 #include "src/snes/impls/tr/tr.h"
26 #include "log.h"
27 #include "opt_update.h"
28 #include "gen_update.h"
29
30 /* ----------------------------------------------------------------------------
31 * QDDM_Solver_L1E: constructor
32 */
33 QDDM_Solver_L1E::QDDM_Solver_L1E():N(0)
34 {
35 its = 0;
36 norm = 0;
37 mf_flg = PETSC_FALSE;
38 }
39
40
41 /* ----------------------------------------------------------------------------
42 * form_function_pn: This function setup DDM equation F(x)=0
43 */
44 void QDDM_Solver_L1E::form_function_pn_1Q(PetscScalar *x,PetscScalar *f)
45 {
46 // compute flux along triangle edges. semiconductor zone only
47 for(int z=0;z<zone_num;z++)
48 if(zonedata[z]->material_type==Semiconductor)
49 {
50 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
51 for(int i=0;i<zone[z].datri.size();i++)
52 {
53 Tri *ptri = &zone[z].datri[i];
54 pzonedata->F1Q_Tri_qddm(ptri,x,f,zofs);
55 }
56 }
57
58 for(int z=0;z<zone_num;z++)
59 {
60 //semiconductor zone
61 if(zonedata[z]->material_type==Semiconductor)
62 {
63 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
64
65 // process electrode.
66 for(int i=0;i<pzonedata->electrode.size();i++)
67 {
68 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
69 pzonedata->F1Q_om_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
70 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
71 pzonedata->F1Q_stk_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
72 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
73 pzonedata->F1Q_ins_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
74 }
75
76 // process cell variables and boundaries.
77 for(int i=0;i<zone[z].davcell.size();i++)
78 {
79 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
80 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
81 pzonedata->F1Q_qddm_inner(i,x,f,ODE_formula,zofs);
82 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
83 pzonedata->F1Q_qddm_ombc(i,x,f,ODE_formula,zofs,bc);
84 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
85 pzonedata->F1Q_qddm_stkbc(i,x,f,ODE_formula,zofs,bc);
86 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
87 pzonedata->F1Q_qddm_insulator_gate(i,x,f,ODE_formula,zofs,bc);
88 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
89 {
90 InsulatorInterfaceBC *pbc;
91 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
92 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
93 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
94 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
95 pzonedata->F1Q_qddm_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
96 }
97 else //prevent some mistakes caused by inexact bc
98 pzonedata->F1Q_qddm_inner(i,x,f,ODE_formula,zofs);
99 }
100 }
101
102 // Insulator zone
103 if(zonedata[z]->material_type==Insulator)
104 {
105 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
106 pzonedata->F1Q_efield_update(x,zofs,bc,zonedata);
107 for(int i=0;i<zone[z].davcell.size();i++)
108 {
109 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
110 if(!pcell->bc_index || bc[pcell->bc_index-1].BCType==NeumannBoundary)
111 pzonedata->F1Q_ddm_inner(i,x,f,ODE_formula,zofs);
112 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==GateContact)
113 pzonedata->F1Q_ddm_gatebc(i,x,f,ODE_formula,zofs,bc);
114 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
115 pzonedata->F1Q_ddm_chargebc(i,x,f,ODE_formula,zofs,bc);
116 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
117 {
118 ElectrodeInsulatorBC *pbc;
119 pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
120 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
121 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
122 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
123 pzonedata->F1Q_ddm_electrode_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
124 }
125 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==InsulatorInterface)
126 {
127 InsulatorInterfaceBC *pbc;
128 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
129 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
130 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
131 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
132 pzonedata->F1Q_ddm_semiconductor_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
133 }
134 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
135 {
136 InsulatorInsulatorBC *pbc;
137 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
138 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
139 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
140 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
141 pzonedata->F1Q_ddm_insulator_insulator_interface(i,x,f,ODE_formula,zofs,bc,pz,n_node);
142 }
143 else //prevent some mistakes caused by inexact bc
144 pzonedata->F1Q_ddm_inner(i,x,f,ODE_formula,zofs);
145 }
146 for(int i=0;i<pzonedata->electrode.size();i++)
147 {
148 if(bc[pzonedata->electrode[i]].BCType==GateContact)
149 pzonedata->F1Q_gate_electrode(i,x,f,ODE_formula,zofs,bc,DeviceDepth);
150 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
151 pzonedata->F1Q_charge_electrode(i,x,f,ODE_formula,zofs,bc);
152 }
153 }
154
155 // Electrode zone
156 if(zonedata[z]->material_type==Conductor)
157 {
158 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
159 for(int i=0;i<zone[z].davcell.size();i++)
160 {
161 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
162 if(!pcell->bc_index || bc[pcell->bc_index-1].BCType==NeumannBoundary)
163 pzonedata->F1Q_ddm_inner(i,x,f,ODE_formula,zofs);
164 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
165 pzonedata->F1Q_ddm_inner(i,x,f,ODE_formula,zofs);
166 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==OhmicContact)
167 {
168 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
169 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
170 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
171 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
172 pzonedata->F1Q_ddm_om_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
173 }
174 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==SchottkyContact)
175 {
176 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
177 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
178 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
179 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
180 pzonedata->F1Q_ddm_om_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
181 }
182 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==GateContact)
183 {
184 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
185 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
186 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
187 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
188 pzonedata->F1Q_ddm_gate_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
189 }
190 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
191 {
192 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
193 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
194 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
195 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
196 pzonedata->F1Q_ddm_charge_contact(i,x,f,ODE_formula,zofs,bc,pz,n_node);
197 }
198 else //prevent some mistakes caused by inexact bc
199 pzonedata->F1Q_ddm_inner(i,x,f,ODE_formula,zofs);
200 }
201 }
202
203 }
204 }
205
206
207 /* ----------------------------------------------------------------------------
208 * form_jacobian_pn: This function setup jacobian matrix F'(x) of DDM equation F(x)=0
209 * dual carrier edition
210 */
211 void QDDM_Solver_L1E::form_jacobian_pn_1Q(PetscScalar *x, Mat *jac, Mat *jtmp)
212 {
213
214 //Compute Jacobian entries for flux along triangle edges.
215 for(int z=0;z<zone_num;z++)
216 {
217 if(zonedata[z]->material_type==Semiconductor)
218 {
219 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
220
221 // compute flux Jacobian along triangle edges
222 for(int i=0;i<zone[z].datri.size();i++)
223 {
224 Tri *ptri = &zone[z].datri[i];
225 pzonedata->J1Q_Tri_qddm(ptri,x,jtmp,zofs);
226 }
227 }
228 }
229 MatAssemblyBegin(*jtmp,MAT_FINAL_ASSEMBLY);
230 MatAssemblyEnd(*jtmp,MAT_FINAL_ASSEMBLY);
231
232 //Compute Jacobian entries and insert into matrix.
233 for(int z=0;z<zone_num;z++)
234 {
235 if(zonedata[z]->material_type==Semiconductor)
236 {
237 SMCZone *pzonedata= dynamic_cast< SMCZone * >(zonedata[z]);
238
239 // process electrode.
240 for(int i=0;i<pzonedata->electrode.size();i++)
241 {
242 if(bc[pzonedata->electrode[i]].BCType==OhmicContact)
243 pzonedata->J1Q_om_electrode(i,x,jac,&JTmp,ODE_formula,zofs,bc,DeviceDepth);
244 if(bc[pzonedata->electrode[i]].BCType==SchottkyContact)
245 pzonedata->J1Q_stk_electrode(i,x,jac,&JTmp,ODE_formula,zofs,bc,DeviceDepth);
246 if(bc[pzonedata->electrode[i]].BCType==InsulatorContact)
247 pzonedata->J1Q_ins_electrode(i,x,jac,&JTmp,ODE_formula,zofs,bc,DeviceDepth);
248 }
249
250 // process cell variables and boundaries.
251 for(int i=0;i<zone[z].davcell.size();i++)
252 {
253 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
254 SemiData *pfs = &pzonedata->fs[i];
255 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
256 pzonedata->J1Q_qddm_inner(i,x,jac,&JTmp,ODE_formula,zofs);
257 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==OhmicContact)
258 pzonedata->J1Q_qddm_ombc(i,x,jac,&JTmp,ODE_formula,zofs,bc);
259 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==SchottkyContact)
260 pzonedata->J1Q_qddm_stkbc(i,x,jac,&JTmp,ODE_formula,zofs,bc);
261 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorContact)
262 pzonedata->J1Q_qddm_insulator_gate(i,x,jac,&JTmp,ODE_formula,zofs,bc);
263 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
264 {
265 InsulatorInterfaceBC *pbc;
266 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
267 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
268 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
269 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
270 pzonedata->J1Q_qddm_interface(i,x,jac,&JTmp,ODE_formula,zofs,bc,pz,n_node);
271 }
272 else //prevent some mistakes caused by inexact bc
273 pzonedata->J1Q_qddm_inner(i,x,jac,&JTmp,ODE_formula,zofs);
274 }
275 }
276
277 // Insulator zone
278 if(zonedata[z]->material_type==Insulator)
279 {
280 ISZone *pzonedata= dynamic_cast< ISZone * >(zonedata[z]);
281
282 for(int i=0;i<zone[z].davcell.size();i++)
283 {
284
285 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
286 if(!pcell->bc_index||bc[pcell->bc_index-1].BCType==NeumannBoundary)
287 pzonedata->J1Q_ddm_inner(i,x,jac,ODE_formula,zofs);
288 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==GateContact)
289 pzonedata->J1Q_ddm_gatebc(i,x,jac,ODE_formula,zofs,bc);
290 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==ChargedContact)
291 pzonedata->J1Q_ddm_chargebc(i,x,jac,ODE_formula,zofs,bc);
292 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
293 {
294 ElectrodeInsulatorBC *pbc;
295 pbc = dynamic_cast<ElectrodeInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
296 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
297 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
298 ElZone * pz = dynamic_cast<ElZone *>(zonedata[n_zone]);
299 pzonedata->J1Q_ddm_electrode_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
300 }
301 else if(pcell->bc_index&&bc[pcell->bc_index-1].BCType==InsulatorInterface)
302 {
303 InsulatorInterfaceBC *pbc;
304 pbc = dynamic_cast<InsulatorInterfaceBC*>(bc.Get_pointer(pcell->bc_index-1));
305 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
306 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
307 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
308 pzonedata->J1Q_ddm_semiconductor_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
309 }
310 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Insulator_Insulator)
311 {
312 InsulatorInsulatorBC *pbc;
313 pbc = dynamic_cast<InsulatorInsulatorBC*>(bc.Get_pointer(pcell->bc_index-1));
314 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
315 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
316 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
317 pzonedata->J1Q_ddm_insulator_insulator_interface(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
318 }
319 else //prevent some mistakes caused by inexact bc
320 pzonedata->J1Q_ddm_inner(i,x,jac,ODE_formula,zofs);
321 }
322 for(int i=0;i<pzonedata->electrode.size();i++)
323 {
324 if(bc[pzonedata->electrode[i]].BCType==GateContact)
325 pzonedata->J1Q_gate_electrode(i,x,jac,ODE_formula,zofs,bc,DeviceDepth);
326 if(bc[pzonedata->electrode[i]].BCType==ChargedContact)
327 pzonedata->J1Q_charge_electrode(i,x,jac,ODE_formula,zofs,bc);
328 }
329 }
330
331 // Electrode zone
332 if(zonedata[z]->material_type==Conductor)
333 {
334 ElZone *pzonedata= dynamic_cast< ElZone * >(zonedata[z]);
335 for(int i=0;i<zone[z].davcell.size();i++)
336 {
337 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
338 if(!pcell->bc_index || bc[pcell->bc_index-1].BCType==NeumannBoundary)
339 pzonedata->J1Q_ddm_inner(i,x,jac,ODE_formula,zofs);
340 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==IF_Electrode_Insulator)
341 pzonedata->J1Q_ddm_inner(i,x,jac,ODE_formula,zofs);
342 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==OhmicContact)
343 {
344 OhmicBC *pbc = dynamic_cast<OhmicBC*>(bc.Get_pointer(pcell->bc_index-1));
345 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
346 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
347 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
348 pzonedata->J1Q_ddm_om_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
349 }
350 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==SchottkyContact)
351 {
352 SchottkyBC *pbc = dynamic_cast<SchottkyBC*>(bc.Get_pointer(pcell->bc_index-1));
353 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
354 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
355 SMCZone * pz = dynamic_cast<SMCZone *>(zonedata[n_zone]);
356 pzonedata->J1Q_ddm_stk_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
357 }
358 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==GateContact)
359 {
360 GateBC *pbc = dynamic_cast<GateBC*>(bc.Get_pointer(pcell->bc_index-1));
361 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
362 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
363 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
364 pzonedata->J1Q_ddm_gate_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
365 }
366 else if(pcell->bc_index && bc[pcell->bc_index-1].BCType==ChargedContact)
367 {
368 ChargedContactBC *pbc = dynamic_cast<ChargedContactBC*>(bc.Get_pointer(pcell->bc_index-1));
369 int n_zone = pbc->pinterface->Find_neighbor_zone_index(z);
370 int n_node = pbc->pinterface->Find_neighbor_node_index(z,i);
371 ISZone * pz = dynamic_cast<ISZone *>(zonedata[n_zone]);
372 pzonedata->J1Q_ddm_charge_contact(i,x,jac,ODE_formula,zofs,bc,pz,n_node);
373 }
374 else //prevent some mistakes caused by inexact bc
375 pzonedata->J1Q_ddm_inner(i,x,jac,ODE_formula,zofs);
376 }
377 }
378
379 }
380 }
381
382
383
384
385 /* ----------------------------------------------------------------------------
386 * SNES_form_function_pn_1: wrap function for petsc nonlinear solver
387 */
388 inline PetscErrorCode SNES_form_function_pn_1Q(SNES snes, Vec x,Vec f,void *dummy)
389 {
390 PetscScalar *xx,*ff;
391 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)dummy;
392 VecZeroEntries(f);
393 //Get pointers to vector data.
394 VecGetArray(x,&xx);
395 VecGetArray(f,&ff);
396
397 ps->form_function_pn_1Q(xx,ff);
398 ps->error_norm_pn_1Q(xx,ff);
399
400 //Restore vectors
401 VecRestoreArray(x,&xx);
402 VecRestoreArray(f,&ff);
403 VecNorm(f,NORM_2,&ps->norm);
404 return 0;
405 }
406
407
408 /* ----------------------------------------------------------------------------
409 * SNES_form_jacobian_pn_1: wrap function for petsc nonlinear solver
410 */
411 inline PetscErrorCode SNES_form_jacobian_pn_1Q(SNES snes, Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
412 {
413 PetscScalar *xx;
414 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)dummy;
415 //Get pointer to vector data
416 VecGetArray(x,&xx);
417 //clear old matrix
418 MatZeroEntries(*jac);
419 MatZeroEntries(ps->JTmp);
420 //build matrix here
421 ps->form_jacobian_pn_1Q(xx,jac,&ps->JTmp);
422
423 *flag = SAME_NONZERO_PATTERN;
424 //Restore vector
425 VecRestoreArray(x,&xx);
426 //Assemble matrix
427 MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
428 MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
429 //MatView(*jac,PETSC_VIEWER_DRAW_WORLD);
430 //getchar();
431 return 0;
432 }
433
434 /* ----------------------------------------------------------------------------
435 * SNESMF_form_jacobian_pn_1: wrap function for petsc nonlinear solver with
436 * matrix-free method
437 */
438 inline PetscErrorCode SNESMF_form_jacobian_pn_1Q(SNES snes, Vec x,Mat *A,Mat *B,MatStructure *flag,void *dummy)
439 {
440 PetscScalar *xx;
441 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)dummy;
442 //Get pointer to vector data
443 VecGetArray(x,&xx);
444 //clear old matrix
445 MatZeroEntries(*B);
446 MatZeroEntries(ps->JTmp);
447 //build matrix here
448 ps->form_jacobian_pn_1Q(xx,B,&ps->JTmp);
449 *flag = SAME_NONZERO_PATTERN;
450 //Restore vector
451 VecRestoreArray(x,&xx);
452 //Assemble matrix
453 MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
454 MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
455 MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
456 MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
457 //MatView(*B,PETSC_VIEWER_DRAW_WORLD);
458 return 0;
459 }
460
461
462 /* ----------------------------------------------------------------------------
463 * error_norm_pn_1Q: This function compute X and RHS error norms.
464 */
465 void QDDM_Solver_L1E:: error_norm_pn_1Q(PetscScalar *x,PetscScalar *f)
466 {
467 // do clear
468 possion_norm = 0;
469 elec_continuty_norm = 0;
470 hole_continuty_norm = 0;
471 elec_quantum_norm = 0;
472 hole_quantum_norm = 0;
473 electrode_norm = 0;
474
475 int offset=0;
476 for(int z=0;z<zone_num;z++)
477 {
478 if(zonedata[z]->material_type==Semiconductor)
479 {
480 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
481 for(int i=0;i<zone[z].davcell.size();i++)
482 {
483 possion_norm += f[offset+0]*f[offset+0];
484 elec_continuty_norm += f[offset+1]*f[offset+1];
485 hole_continuty_norm += f[offset+2]*f[offset+2];
486 elec_quantum_norm += f[offset+3]*f[offset+3];
487 hole_quantum_norm += f[offset+4]*f[offset+4];
488 offset += 5;
489 }
490 for(int i=0;i<pzonedata->electrode.size();i++)
491 {
492 electrode_norm += f[offset]*f[offset];
493 offset += 1;
494 }
495 }
496 if(zonedata[z]->material_type==Insulator)
497 {
498 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
499 for(int i=0;i<zone[z].davcell.size();i++)
500 {
501 possion_norm += f[offset]*f[offset];
502 offset += 1;
503 }
504 for(int i=0;i<pzonedata->electrode.size();i++)
505 {
506 electrode_norm += f[offset]*f[offset];
507 offset += 1;
508 }
509 }
510 if(zonedata[z]->material_type==Conductor)
511 {
512 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
513 for(int i=0;i<zone[z].davcell.size();i++)
514 {
515 possion_norm += f[offset]*f[offset];
516 offset += 1;
517 }
518 }
519 }
520
521 possion_norm = sqrt(possion_norm);
522 elec_continuty_norm = sqrt(elec_continuty_norm);
523 hole_continuty_norm = sqrt(hole_continuty_norm);
524 elec_quantum_norm = sqrt(elec_quantum_norm);
525 hole_quantum_norm = sqrt(hole_quantum_norm);
526 electrode_norm = sqrt(electrode_norm);
527
528 }
529
530
531 inline PetscErrorCode KSPConvergenceTest_L1Q(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason* reason,void *dummy)
532 {
533 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)dummy;
534 ksp->rtol = 1e-10*ps->N;
535 ksp->abstol= PetscMax(1e-6*ps->norm,1e-15*ps->N);
536 return KSPDefaultConverged(ksp,n,rnorm,reason,0);
537 }
538
539
540 inline PetscErrorCode LineSearch_ConvergenceTest_L1Q(SNES snes,PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm,
541 SNESConvergedReason *reason,void *dummy)
542 {
543 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)dummy;
544
545 *reason = SNES_CONVERGED_ITERATING;
546 if (!it)
547 {
548 snes->ttol = fnorm*snes->rtol;
549 gss_log.string_buf()<<" "<<"its\t"<<"| Eq(V) | "<<"| Eq(n) | "<<"| Eq(p) | "<<"|Eq(Qn)| "<<"|Eq(Qp)| "<<"|delta x|\n"
550 <<"-----------------------------------------------------------------------------\n";
551 gss_log.record();
552 }
553 gss_log.string_buf().precision(2);
554 gss_log.string_buf()<<" "<<it<<"\t"
555 <<ps->possion_norm<<" "
556 <<ps->elec_continuty_norm<<" "
557 <<ps->hole_continuty_norm<<" "
558 <<ps->elec_quantum_norm<<" "
559 <<ps->hole_quantum_norm<<" "
560 <<pnorm<<"\n";
561 gss_log.record();
562 gss_log.string_buf().precision(6);
563
564 if (fnorm != fnorm)
565 {
566 *reason = SNES_DIVERGED_FNORM_NAN;
567 }
568 else if (ps->possion_norm < ps->possion_abs_toler &&
569 ps->elec_continuty_norm < ps->elec_continuty_abs_toler &&
570 ps->hole_continuty_norm < ps->hole_continuty_abs_toler &&
571 ps->elec_quantum_norm < ps->elec_quantum_abs_toler &&
572 ps->hole_quantum_norm < ps->hole_quantum_abs_toler &&
573 ps->electrode_norm < ps->electrode_abs_toler)
574 {
575 *reason = SNES_CONVERGED_FNORM_ABS;
576 }
577 else if (snes->nfuncs >= snes->max_funcs)
578 {
579 *reason = SNES_DIVERGED_FUNCTION_COUNT;
580 }
581
582 if (it && !*reason)
583 {
584 if (fnorm <= snes->ttol)
585 {
586 *reason = SNES_CONVERGED_FNORM_RELATIVE;
587 }
588 else if (pnorm < ps->relative_toler &&
589 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
590 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
591 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
592 ps->elec_quantum_norm < ps->toler_relax*ps->elec_quantum_abs_toler &&
593 ps->hole_quantum_norm < ps->toler_relax*ps->hole_quantum_abs_toler &&
594 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler)
595 {
596 *reason = SNES_CONVERGED_PNORM_RELATIVE;
597 }
598 }
599
600 if(*reason>0)
601 {
602 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
603 <<" "<<SNESConvergedReasons[*reason]<<" *residual norm = "<<fnorm<<"\n\n\n";
604 gss_log.record();
605 }
606
607 return(0);
608 }
609
610
611 inline PetscErrorCode TrustRegion_ConvergenceTest_L1Q(SNES snes,PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm,
612 SNESConvergedReason *reason,void *dummy)
613 {
614 SNES_TR *neP = (SNES_TR *)snes->data;
615 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)dummy;
616
617 if (!ps->its)
618 {
619 gss_log.string_buf()<<" "<<"its\t"<<"| Eq(V) | "<<"| Eq(n) | "<<"| Eq(p) | "<<"|Eq(Qn)| "<<"|Eq(Qp)| "<<"|delta x|\n"
620 <<"-----------------------------------------------------------------------------\n";
621 gss_log.record();
622 }
623 gss_log.string_buf().precision(2);
624 gss_log.string_buf()<<" "<<it<<"\t"
625 <<ps->possion_norm<<" "
626 <<ps->elec_continuty_norm<<" "
627 <<ps->hole_continuty_norm<<" "
628 <<ps->elec_quantum_norm<<" "
629 <<ps->hole_quantum_norm<<" "
630 <<pnorm<<"\n";
631 gss_log.record();
632 gss_log.string_buf().precision(6);
633
634 *reason = SNES_CONVERGED_ITERATING;
635
636 if (!it)
637 {
638 snes->ttol = fnorm*snes->rtol;
639 }
640
641 if (fnorm != fnorm)
642 {
643 *reason = SNES_DIVERGED_FNORM_NAN;
644 }
645 else if (neP->delta < xnorm * snes->deltatol)
646 {
647 if( ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
648 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
649 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
650 ps->elec_quantum_norm < ps->toler_relax*ps->elec_quantum_abs_toler &&
651 ps->hole_quantum_norm < ps->toler_relax*ps->hole_quantum_abs_toler &&
652 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler)
653 {
654 *reason = SNES_CONVERGED_TR_DELTA;
655 }
656 else
657 {
658 *reason = SNES_DIVERGED_LOCAL_MIN;
659 }
660 }
661 else if (ps->possion_norm < ps->possion_abs_toler &&
662 ps->elec_continuty_norm < ps->elec_continuty_abs_toler &&
663 ps->hole_continuty_norm < ps->hole_continuty_abs_toler &&
664 ps->elec_quantum_norm < ps->elec_quantum_abs_toler &&
665 ps->hole_quantum_norm < ps->hole_quantum_abs_toler &&
666 ps->electrode_norm < ps->electrode_abs_toler)
667 {
668 *reason = SNES_CONVERGED_FNORM_ABS;
669 }
670 else if (snes->nfuncs >= snes->max_funcs)
671 {
672 *reason = SNES_DIVERGED_FUNCTION_COUNT;
673 }
674
675 if (it && !*reason)
676 {
677 if (fnorm <= snes->ttol)
678 {
679 *reason = SNES_CONVERGED_FNORM_RELATIVE;
680 }
681 else if (pnorm < ps->relative_toler &&
682 ps->possion_norm < ps->toler_relax*ps->possion_abs_toler &&
683 ps->elec_continuty_norm < ps->toler_relax*ps->elec_continuty_abs_toler &&
684 ps->hole_continuty_norm < ps->toler_relax*ps->hole_continuty_abs_toler &&
685 ps->elec_quantum_norm < ps->toler_relax*ps->elec_quantum_abs_toler &&
686 ps->hole_quantum_norm < ps->toler_relax*ps->hole_quantum_abs_toler &&
687 ps->electrode_norm < ps->toler_relax*ps->electrode_abs_toler)
688 {
689 *reason = SNES_CONVERGED_PNORM_RELATIVE;
690 }
691 }
692
693 if (snes->nfuncs >= snes->max_funcs)
694 {
695 *reason = SNES_DIVERGED_FUNCTION_COUNT;
696 }
697 if(*reason>0)
698 {
699 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
700 <<" "<<SNESConvergedReasons[*reason]<<" *residual norm = "<<fnorm<<"\n\n\n";
701 gss_log.record();
702 }
703
704 return(0);
705 }
706
707 PetscErrorCode LimitorByBankRose_L1Q(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
708 {
709 int it;
710 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)checkctx;
711 PetscScalar *ww;
712 VecGetArray(w,&ww);
713
714 SNESGetIterationNumber(snes,&it);
715 static PetscScalar K;
716 if(it==0) K=0;
717
718 Vec gk0,gk1;
719 int flag=0;
720 VecDuplicate(x,&gk0);
721 VecDuplicate(x,&gk1);
722 SNESComputeFunction(snes,x,gk0);
723 PetscScalar gk0_norm;
724 VecNorm(gk0,NORM_2,&gk0_norm);
725 while(1)
726 {
727 PetscScalar tk=1.0/(1+K*gk0_norm);
728 if(tk<1e-4) break;
729 //w=x-tk*y;
730 flag=1;
731 VecWAXPY(w,-tk,y,x);
732 SNESComputeFunction(snes,w,gk1);
733 PetscScalar gk1_norm;
734 VecNorm(gk1,NORM_2,&gk1_norm);
735 if((1-gk1_norm/gk0_norm)<0.9*tk)
736 {
737 if(K==0) K=1;
738 else K*=10;
739 }
740 else
741 {
742 K/=10;
743 break;
744 }
745 }
746
747 //prevent negative carrier density
748 int offset=0;
749 for(int z=0;z<ps->zone_num;z++)
750 {
751 if(ps->zonedata[z]->material_type==Semiconductor)
752 {
753 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
754 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
755 {
756 if (ww[offset+1] <0 ) {ww[offset+1]=fabs(ww[offset+1]);flag=1;}
757 if (ww[offset+2] <0 ) {ww[offset+2]=fabs(ww[offset+2]);flag=1;}
758 offset += 5;
759 }
760 offset += pzonedata->electrode.size();
761 }
762 else if(ps->zonedata[z]->material_type==Insulator)
763 {
764 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
765 offset += pzonedata->pzone->davcell.size();
766 offset += pzonedata->electrode.size();
767 }
768 else if(ps->zonedata[z]->material_type==Conductor)
769 {
770 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
771 offset += pzonedata->pzone->davcell.size();
772 }
773 }
774
775 if(flag)
776 {
777 *changed_y = PETSC_FALSE;
778 *changed_w = PETSC_TRUE;
779 }
780 VecDestroy(gk0);
781 VecDestroy(gk1);
782 VecRestoreArray(w,&ww);
783 return(0);
784 }
785
786
787 PetscErrorCode LimitorByPotential_L1Q(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
788 {
789 PetscScalar *xx;
790 PetscScalar *yy;
791 PetscScalar *ww;
792 VecGetArray(x,&xx);
793 VecGetArray(y,&yy);
794 VecGetArray(w,&ww);
795 PetscScalar dV_max=0;
796
797 //search for dV_max;
798 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)checkctx;
799 int offset=0;
800 for(int z=0;z<ps->zone_num;z++)
801 {
802 if(ps->zonedata[z]->material_type==Semiconductor)
803 {
804 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
805 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
806 {
807 if(fabs(yy[offset])>dV_max) {dV_max=fabs(yy[offset]); }
808 if(fabs(yy[offset+3])>dV_max) {dV_max=fabs(yy[offset+3]); }
809 if(fabs(yy[offset+4])>dV_max) {dV_max=fabs(yy[offset+4]); }
810 offset += 5;
811 }
812 offset += pzonedata->electrode.size();
813 }
814 else if(ps->zonedata[z]->material_type==Insulator)
815 {
816 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
817 offset += pzonedata->pzone->davcell.size();
818 offset += pzonedata->electrode.size();
819 }
820 else if(ps->zonedata[z]->material_type==Conductor)
821 {
822 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
823 offset += pzonedata->pzone->davcell.size();
824 }
825 }
826 if(dV_max<1e-8) dV_max=1e-8;
827
828 //do logarithmic potential damping;
829 PetscScalar Vt=0.026*ps->LatticeTemp;
830 PetscScalar f = log(1+dV_max/Vt)/(dV_max/Vt);
831
832 //prevent negative carrier density
833 offset=0;
834 for(int z=0;z<ps->zone_num;z++)
835 {
836 if(ps->zonedata[z]->material_type==Semiconductor)
837 {
838 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
839 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
840 {
841 ww[offset] = xx[offset] - f*yy[offset];
842 if (ww[offset+1]<0.0) ww[offset+1]=fabs(ww[offset+1]);
843 if (ww[offset+2]<0.0) ww[offset+2]=fabs(ww[offset+2]);
844 ww[offset+3] = xx[offset+3] - f*yy[offset+3];
845 ww[offset+4] = xx[offset+4] - f*yy[offset+4];
846 offset += 5;
847 }
848 offset += pzonedata->electrode.size();
849 }
850 else if(ps->zonedata[z]->material_type==Insulator)
851 {
852 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
853 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
854 {
855 ww[offset] = xx[offset] - f*yy[offset];
856 offset += 1;
857 }
858 offset += pzonedata->electrode.size();
859 }
860 else if(ps->zonedata[z]->material_type==Conductor)
861 {
862 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
863 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
864 {
865 ww[offset] = xx[offset] - f*yy[offset];
866 offset += 1;
867 }
868 }
869 }
870
871 VecRestoreArray(x,&xx);
872 VecRestoreArray(y,&yy);
873 VecRestoreArray(w,&ww);
874
875 *changed_y = PETSC_FALSE;
876 *changed_w = PETSC_TRUE;
877 return(0);
878 }
879
880
881
882
883 PetscErrorCode LimitorNonNegativeCarrier_L1Q(SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
884 {
885 PetscScalar *xx;
886 PetscScalar *yy;
887 PetscScalar *ww;
888 VecGetArray(x,&xx);
889 VecGetArray(y,&yy);
890 VecGetArray(w,&ww);
891 PetscScalar WorstRatio=0.0;
892
893 QDDM_Solver_L1E *ps = (QDDM_Solver_L1E*)checkctx;
894 int offset=0;
895 for(int z=0;z<ps->zone_num;z++)
896 {
897 if(ps->zonedata[z]->material_type==Semiconductor)
898 {
899 SMCZone *pzonedata = dynamic_cast< SMCZone * >(ps->zonedata[z]);
900 for(int i=0;i<pzonedata->pzone->davcell.size();i++)
901 {
902 PetscScalar r;
903 if (ww[offset+1]<0.0) ww[offset+1]=fabs(ww[offset+1]);
904 if (ww[offset+2]<0.0) ww[offset+2]=fabs(ww[offset+2]);
905 offset += 5;
906 }
907 offset += pzonedata->electrode.size();
908 }
909 else if(ps->zonedata[z]->material_type==Insulator)
910 {
911 ISZone *pzonedata = dynamic_cast< ISZone * >(ps->zonedata[z]);
912 offset += pzonedata->pzone->davcell.size();
913 offset += pzonedata->electrode.size();
914 }
915 else if(ps->zonedata[z]->material_type==Conductor)
916 {
917 ElZone *pzonedata = dynamic_cast< ElZone * >(ps->zonedata[z]);
918 offset += pzonedata->pzone->davcell.size();
919 }
920 }
921 VecRestoreArray(x,&xx);
922 VecRestoreArray(y,&yy);
923 VecRestoreArray(w,&ww);
924
925 *changed_y = PETSC_FALSE;
926 *changed_w = PETSC_TRUE;
927 return(0);
928 }
929
930
931 /* ----------------------------------------------------------------------------
932 * QDDM_Solver_L1E::init_solver: This function do initial setup for nonlinear solver
933 */
934 int QDDM_Solver_L1E::init_solver(SolveDefine &sv)
935 {
936 gss_log.string_buf()<<"DG-DDM solver Level 1E init...";
937
938 //set Tolerances
939 relative_toler = sv.relative_toler;
940 toler_relax = sv.toler_relax;
941 possion_abs_toler = sv.possion_abs_toler;
942 elec_continuty_abs_toler = sv.elec_continuty_abs_toler;
943 hole_continuty_abs_toler = sv.hole_continuty_abs_toler;
944 electrode_abs_toler = sv.electrode_abs_toler;
945 elec_quantum_abs_toler = sv.elec_quantum_abs_toler;
946 hole_quantum_abs_toler = sv.hole_quantum_abs_toler;
947
948 //if user defined a matrix-free method
949 PetscOptionsHasName(PETSC_NULL,"-snes_mf_operator",&mf_flg);
950
951 // compute the scale of problem
952 zofs.resize(zone_num+1);
953 zofs[0] = 0;
954 for(int i=0;i<zone_num;i++)
955 {
956 if(zonedata[i]->material_type==Semiconductor) //semiconductor zone
957 {
958 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
959 N += 5*pzonedata->pzone->davcell.size();
960 N += pzonedata->electrode.size(); //additional equation for electrode
961 zofs[i+1] = N;
962 }
963 else if(zonedata[i]->material_type==Insulator) //Insulator zone
964 {
965 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
966 N += zone[i].davcell.size();
967 N += pzonedata->electrode.size(); //additional equation for electrode
968 zofs[i+1] = N;
969 }
970 else if(zonedata[i]->material_type==Conductor) //Electrode zone
971 {
972 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
973 N += zone[i].davcell.size();
974 zofs[i+1] = N;
975 }
976 else
977 {
978 zofs[i+1] = N;
979 }
980 }
981
982 VecCreateSeq(PETSC_COMM_SELF,N,&x);
983 VecDuplicate(x,&r);
984 //Evaluate initial guess,importnat for newton solver
985 PetscScalar init_value[5];
986 int index[5];
987 int offset=0;
988 for(int z=0;z<zone_num;z++)
989 {
990 if(zonedata[z]->material_type==Semiconductor)
991 {
992 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
993 for(int i=0;i<zone[z].davcell.size();i++)
994 {
995 index[0] = offset+0;
996 index[1] = offset+1;
997 index[2] = offset+2;
998 index[3] = offset+3;
999 index[4] = offset+4;
1000 init_value[0] = pzonedata->fs[i].P;
1001 init_value[1] = pzonedata->fs[i].n;
1002 init_value[2] = pzonedata->fs[i].p;
1003 init_value[3] = pzonedata->fs[i].Eqc;
1004 init_value[4] = pzonedata->fs[i].Eqv;
1005 VecSetValues(x,5,index,init_value,INSERT_VALUES);
1006 offset += 5;
1007 }
1008 for(int j=0;j<pzonedata->electrode.size();j++)
1009 {
1010 VecSetValue(x,offset,bc.Get_pointer(pzonedata->electrode[j])->Get_Potential(),INSERT_VALUES);
1011 offset += 1;
1012 }
1013 }
1014
1015 if(zonedata[z]->material_type==Insulator)
1016 {
1017 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1018 for(int i=0;i<zone[z].davcell.size();i++)
1019 {
1020 VecSetValue(x,offset,pzonedata->fs[i].P,INSERT_VALUES);
1021 offset += 1;
1022 }
1023 for(int j=0;j<pzonedata->electrode.size();j++)
1024 {
1025 VecSetValue(x,offset,bc.Get_pointer(pzonedata->electrode[j])->Get_Potential(),INSERT_VALUES);
1026 offset += 1;
1027 }
1028 }
1029
1030 if(zonedata[z]->material_type==Conductor)
1031 {
1032 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1033 for(int i=0;i<zone[z].davcell.size();i++)
1034 {
1035 VecSetValue(x,offset,pzonedata->fs[i].P,INSERT_VALUES);
1036 offset += 1;
1037 }
1038 }
1039
1040 }
1041 //
1042 VecAssemblyBegin(x);
1043 VecAssemblyEnd(x);
1044
1045 // Create Jacobian matrix data structure
1046 // pre-alloc approximate memory
1047 int *nnz = new int[N];
1048 int *p = nnz;
1049 for(int i=0;i<zone_num;i++)
1050 {
1051 //semiconductor zone
1052 if(zonedata[i]->material_type==Semiconductor)
1053 {
1054 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1055 for(int j=0;j<pzonedata->pzone->davcell.size();j++)
1056 {
1057 const VoronoiCell* pcell = pzonedata->pzone->davcell.GetPointer(j);
1058 //inner node, alloc exact memory.
1059 if(!pcell->bc_index || bc[pcell->bc_index-1].psegment->interface==-1)
1060 {
1061 *p++ = 5*(pzonedata->pzone->davcell[j].nb_num+1);
1062 *p++ = 5*(pzonedata->pzone->davcell[j].nb_num+1);
1063 *p++ = 5*(pzonedata->pzone->davcell[j].nb_num+1);
1064 *p++ = 5*(pzonedata->pzone->davcell[j].nb_num+1);
1065 *p++ = 5*(pzonedata->pzone->davcell[j].nb_num+1);
1066 }
1067 else //interface node, slightly more than needed.
1068 {
1069 *p++ = 5*(2*pzonedata->pzone->davcell[j].nb_num+4);
1070 *p++ = 5*(2*pzonedata->pzone->davcell[j].nb_num+4);
1071 *p++ = 5*(2*pzonedata->pzone->davcell[j].nb_num+4);
1072 *p++ = 5*(2*pzonedata->pzone->davcell[j].nb_num+4);
1073 *p++ = 5*(2*pzonedata->pzone->davcell[j].nb_num+4);
1074 }
1075 }
1076 for(int j=0;j<pzonedata->electrode.size();j++)
1077 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
1078
1079 }
1080 //Insulator zones, poisson equation only
1081 if(zonedata[i]->material_type==Insulator)
1082 {
1083 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1084 for(int j=0;j<zone[i].davcell.size();j++)
1085 *p++ = 1*(zone[i].davcell[j].nb_num+1);
1086 for(int j=0;j<pzonedata->electrode.size();j++)
1087 *p++ = bc[pzonedata->electrode[j]].psegment->node_array.size()+1;
1088 }
1089 //Electrode zones, poisson equation only
1090 if(zonedata[i]->material_type==Conductor)
1091 {
1092 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[i]);
1093 for(int j=0;j<zone[i].davcell.size();j++)
1094 *p++ = 1*(zone[i].davcell[j].nb_num+1);
1095 }
1096 }
1097 MatCreate(PETSC_COMM_SELF,&J);
1098 MatSetSizes(J,N,N,N,N);
1099 if(sv.LS=="superlu")
1100 {
1101 #ifdef PETSC_HAVE_SUPERLU
1102 MatSetType(J,MATSUPERLU);
1103 #else
1104 MatSetType(J,MATSEQAIJ);
1105 #endif
1106 }
1107 else if(sv.LS=="umfpack")
1108 {
1109 #ifdef PETSC_HAVE_UMFPACK
1110 MatSetType(J,MATUMFPACK);
1111 #else
1112 MatSetType(J,MATSEQAIJ);
1113 #endif
1114 }
1115 else
1116 {
1117 MatSetType(J,MATSEQAIJ);
1118 }
1119 MatSeqAIJSetPreallocation(J,0,nnz);
1120
1121 MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,0,nnz,&JTmp);
1122 if(mf_flg)
1123 MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,0,nnz,&JPrec);
1124
1125 delete [] nnz;
1126
1127
1128 SNESCreate(PETSC_COMM_WORLD,&snes);
1129 SNESSetFunction(snes,r,SNES_form_function_pn_1Q,this);
1130 if(mf_flg)
1131 {
1132 gss_log.string_buf()<<"Matrix-Free...";
1133 SNESSetJacobian(snes,J,JPrec,SNESMF_form_jacobian_pn_1Q,this);
1134 }
1135 else
1136 {
1137 SNESSetJacobian(snes,J,J,SNES_form_jacobian_pn_1Q,this);
1138 }
1139 // set the newton method
1140 SNESSetType(snes,SNESLS); //default method
1141 // the maximum number of iterations
1142 PetscInt maxit = sv.maxit;
1143 if(sv.Type==EQUILIBRIUM) maxit = 10*sv.maxit;
1144
1145 if(sv.NS==LineSearch || sv.NS==Basic)
1146 {
1147 if(sv.NS==LineSearch)
1148 SNESLineSearchSet(snes,SNESLineSearchCubic,PETSC_NULL);
1149 if(sv.NS==Basic)
1150 SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);
1151
1152 if(sv.Damping==DampingBankRose)
1153 SNESLineSearchSetPostCheck(snes,LimitorByBankRose_L1Q,this);
1154 else if(sv.Damping==DampingPotential)
1155 SNESLineSearchSetPostCheck(snes,LimitorByPotential_L1Q,this);
1156 else
1157 SNESLineSearchSetPostCheck(snes,LimitorNonNegativeCarrier_L1Q,this);
1158
1159 SNESSetTolerances(snes,1e-12*N,1e-14,1e-9,maxit,1000);
1160 SNESSetConvergenceTest(snes,LineSearch_ConvergenceTest_L1Q,this);
1161 }
1162 else if(sv.NS==TrustRegion)
1163 {
1164 SNESSetType(snes,SNESTR);
1165 SNESSetTolerances(snes,1e-12*N,1e-14,1e-9,maxit,1000);
1166 SNESSetTrustRegionTolerance(snes,1e-30);
1167 SNESSetConvergenceTest(snes,TrustRegion_ConvergenceTest_L1Q,this);
1168 //must set TR delta0 to a sufficient large value, or TR can't find real solution.
1169 SNES_TR *neP = (SNES_TR*)snes->data;
1170 neP->delta0 = N;
1171 }
1172
1173 SNESGetKSP(snes,&ksp);
1174 KSPGetPC(ksp,&pc);
1175 //KSPSetMonitor(ksp,KSPDefaultMonitor,PETSC_NULL,PETSC_NULL);
1176 if(sv.LS=="lu"||sv.LS=="superlu"||sv.LS=="umfpack")
1177 {
1178 KSPSetType(ksp,KSPPREONLY);
1179 PCSetType(pc,PCLU);
1180 PCFactorSetReuseOrdering(pc,PETSC_TRUE);
1181 PCFactorSetPivoting(pc,1.0);
1182 PCFactorReorderForNonzeroDiagonal(pc,1e-14);
1183 PCFactorSetShiftNonzero(pc,1e-12);
1184 }
1185 else
1186 {
1187 KSPSetType(ksp,sv.LS.c_str());
1188 if(sv.LS=="gmres") KSPGMRESSetRestart(ksp,150);
1189 KSPSetTolerances(ksp,1e-10*N,1e-20*N,PETSC_DEFAULT,N/10);
1190 PCSetType(pc,PCILU);
1191 PCFactorSetLevels(pc,3);
1192 PCFactorSetShiftNonzero(pc,1e-14);
1193 PCFactorReorderForNonzeroDiagonal(pc,1e-12);
1194 }
1195 KSPSetFromOptions(ksp);
1196 SNESSetFromOptions(snes);
1197 gss_log.string_buf()<<"done\n";
1198 gss_log.record();
1199 return 0;
1200 }
1201
1202
1203
1204 /* ----------------------------------------------------------------------------
1205 * QDDM_Solver_L1E::do_solve: This function solve the problem
1206 */
1207 int QDDM_Solver_L1E::do_solve(SolveDefine &sv)
1208 {
1209 if(sv.Type==EQUILIBRIUM)
1210 solve_equ(sv);
1211
1212 else if(sv.Type==STEADYSTATE)
1213 solve_steadystate(sv);
1214
1215 else if(sv.Type==DCSWEEP)
1216 solve_dcsweep(sv);
1217
1218 else if(sv.Type==TRANSIENT)
1219 solve_transient(sv);
1220
1221 else
1222 {
1223 gss_log.string_buf()<<"DG-DDM L1E not support this solver type!\n";
1224 gss_log.record();
1225 }
1226 return 0;
1227
1228 }
1229
1230 /* ----------------------------------------------------------------------------
1231 * QDDM_Solver_L1E::solve_equ: This function compute equilibrium
1232 * all the stimulate source(s) are zero. time step set to inf
1233 */
1234 int QDDM_Solver_L1E::solve_equ(SolveDefine &sv)
1235 {
1236 gss_log.string_buf()<<"Compute equilibrium\n";
1237 gss_log.record();
1238 bc.Clear_Vapp();
1239 bc.Clear_Iapp();
1240 ODE_formula.TimeDependent = false;
1241 ODE_formula.dt = 1e100;
1242 ODE_formula.clock = 0.0;
1243 its = 0;
1244
1245
1246 for(int i=0;i<zone_num;i++)
1247 if(zonedata[i]->material_type == Semiconductor)
1248 {
1249 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1250 pzonedata->HighFieldMobility= false;
1251 pzonedata->ImpactIonization = false;
1252 pzonedata->IIType = sv.IIType;
1253 pzonedata->BandBandTunneling = false;
1254 pzonedata->IncompleteIonization = false;
1255 pzonedata->QuantumMechanical = false;
1256 pzonedata->Fermi = sv.Fermi;
1257 pzonedata->EJModel = sv.EJModel;
1258 pzonedata->QNFactor = sv.QNFactor;
1259 pzonedata->QPFactor = sv.QPFactor;
1260 }
1261 // for equilibrium compute, low field mobility is ok.
1262 solver_pre_compute();
1263
1264 SNESSolve(snes,PETSC_NULL,x);
1265
1266 solution_update();
1267 SNESGetConvergedReason(snes,&reason);
1268 if(reason<0)
1269 {
1270 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
1271 <<" "<<SNESConvergedReasons[reason]<<" !residual norm = "<<norm<<"\n\n\n";
1272 gss_log.record();
1273 }
1274 return 0;
1275 }
1276
1277
1278 /* ----------------------------------------------------------------------------
1279 * QDDM_Solver_L1E::solve_steadystate: This function compute steadystate
1280 * set all the stimulate source(s) at there value of t=0. time step set to inf
1281 */
1282 int QDDM_Solver_L1E::solve_steadystate(SolveDefine &sv)
1283 {
1284
1285 gss_log.string_buf()<<"Compute steady-state\n";
1286 gss_log.record();
1287 bc.Update_Vapp(0);
1288 bc.Update_Iapp(0);
1289 ODE_formula.TimeDependent = false;
1290 ODE_formula.dt = 1e100;
1291 ODE_formula.clock = 0.0;
1292 its = 0;
1293
1294 for(int i=0;i<zone_num;i++)
1295 if(zonedata[i]->material_type == Semiconductor)
1296 {
1297 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1298 pzonedata->HighFieldMobility= sv.HighFieldMobility;
1299 pzonedata->ImpactIonization = sv.ImpactIonization;
1300 pzonedata->IIType = sv.IIType;
1301 pzonedata->BandBandTunneling = sv.BandBandTunneling;
1302 pzonedata->IncompleteIonization = sv.IncompleteIonization;
1303 pzonedata->QuantumMechanical = sv.QuantumMechanical;
1304 pzonedata->Fermi = sv.Fermi;
1305 pzonedata->EJModel = sv.EJModel;
1306 pzonedata->QNFactor = sv.QNFactor;
1307 pzonedata->QPFactor = sv.QPFactor;
1308 }
1309 solver_pre_compute();
1310
1311 SNESSolve(snes,PETSC_NULL,x);
1312 solution_update();
1313 SNESGetConvergedReason(snes,&reason);
1314 if(reason<0)
1315 {
1316 gss_log.string_buf()<<"----------------------------------------------------------------------\n"
1317 <<" "<<SNESConvergedReasons[reason]<<" !residual norm = "<<norm<<"\n\n\n";
1318 gss_log.record();
1319 }
1320 return 0;
1321
1322 }
1323
1324
1325 /* ----------------------------------------------------------------------------
1326 * QDDM_Solver_L1E::solve_dcsweep: This function compute IV curve
1327 * time step set to inf
1328 */
1329 int QDDM_Solver_L1E::solve_dcsweep(SolveDefine &sv)
1330 {
1331 PetscScalar I=0;
1332 FILE *fiv;
1333 PetscScalar current_scale_mA = scale_unit.s_mA;
1334 PetscScalar voltage_scale_V =scale_unit.s_volt;
1335
1336 bc.Update_Vapp(0);
1337 bc.Update_Iapp(0);
1338 optgen_update(0);
1339 ODE_formula.TimeDependent = false;
1340 ODE_formula.dt = 1e100;
1341 ODE_formula.clock = 0.0;
1342
1343 //set which electrodes are required to record IV
1344 for(int i=0;i<zone_num;i++)
1345 {
1346 if(zonedata[i]->material_type == Semiconductor)
1347 {
1348 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1349 pzonedata->HighFieldMobility= sv.HighFieldMobility;
1350 pzonedata->ImpactIonization = sv.ImpactIonization;
1351 pzonedata->IIType = sv.IIType;
1352 pzonedata->BandBandTunneling = sv.BandBandTunneling;
1353 pzonedata->IncompleteIonization = sv.IncompleteIonization;
1354 pzonedata->QuantumMechanical = sv.QuantumMechanical;
1355 pzonedata->Fermi = sv.Fermi;
1356 pzonedata->EJModel = sv.EJModel;
1357 pzonedata->QNFactor = sv.QNFactor;
1358 pzonedata->QPFactor = sv.QPFactor;
1359 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1360 for(int j=0;j<pzonedata->electrode.size();j++)
1361 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1362 {
1363 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1364 sv.Electrode_Record_Index.push_back(k);
1365 }
1366 }
1367 if(zonedata[i]->material_type == Insulator)
1368 {
1369 ISZone * pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1370 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1371 for(int j=0;j<pzonedata->electrode.size();j++)
1372 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1373 {
1374 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1375 sv.Electrode_Record_Index.push_back(k);
1376 }
1377 }
1378 }
1379
1380
1381 // output DC Scan information
1382 if(sv.Electrode_VScan!=-1)
1383 {
1384 if(!sv.Electrode_VScan_Name.size())
1385 {
1386 gss_log.string_buf()<<"No VScan Electrode Specified.\n";
1387 gss_log.record();
1388 return 1;
1389 }
1390 gss_log.string_buf()<<"DC voltage scan from "<<sv.VStart
1391 <<" step "<<sv.VStep
1392 <<" to "<<sv.VStop<<"\n";
1393 gss_log.record();
1394 }
1395 else
1396 {
1397 gss_log.string_buf()<<"DC current scan from "<<sv.IStart/current_scale_mA
1398 <<" step "<<sv.IStep/current_scale_mA
1399 <<" to "<<sv.IStop/current_scale_mA<<"\n";
1400 gss_log.record();
1401 }
1402
1403 // prepare IV file. If no file given, output to screen.
1404 if(!sv.IVFile.empty())
1405 {
1406 if(sv.IVFileAppend)
1407 {
1408 fiv=fopen(sv.IVFile.c_str(),"a");
1409 if(!fiv) fiv=stdout;
1410 }
1411 else
1412 {
1413 fiv=fopen(sv.IVFile.c_str(),"w");
1414 if(!fiv) fiv=stdout;
1415 fprintf(fiv,"#");
1416 for(int j=0;j<sv.Electrode_Record.size();j++)
1417 fprintf(fiv,"Vp(%s)[V] I(%s)[mA] ",
1418 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str(),
1419 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str()
1420 );
1421 fprintf(fiv,"\n");
1422 }
1423 }
1424 else
1425 fiv = stdout;
1426
1427 // begin
1428 if(sv.Electrode_VScan!=-1)
1429 {
1430 PetscScalar V = sv.VStart;
1431 PetscScalar VStep = sv.VStep;
1432 for(int i=0;i<sv.Electrode_VScan_Name.size();i++)
1433 bc.Set_electrode_type(sv.Electrode_VScan_Name[i].c_str(),VoltageBC);
1434 int diverged_retry=0;
1435 probe_open(DCSWEEP_VSCAN);
1436 do
1437 {
1438 its = 0;
1439 gss_log.string_buf()<<"DC Scan: V("<<sv.Electrode_VScan_Name[0];
1440 for(int i=1;i<sv.Electrode_VScan_Name.size();i++)
1441 gss_log.string_buf()<<", "<<sv.Electrode_VScan_Name[i];
1442 gss_log.string_buf()<<") = "<<V/voltage_scale_V<<" V"<<"\n";
1443 gss_log.record();
1444 for(int i=0;i<sv.Electrode_VScan_Name.size();i++)
1445 bc.Set_Vapp(sv.Electrode_VScan_Name[i].c_str(),V);
1446 solver_pre_compute();
1447 SNESSolve(snes,PETSC_NULL,x);
1448 SNESGetConvergedReason(snes,&reason);
1449
1450 if(reason>0 && norm==norm) //ok, converged.
1451 {
1452 diverged_retry=0;
1453 solution_update();
1454 probe(DCSWEEP_VSCAN,V);
1455 if(fabs(VStep) < fabs(sv.VStep))
1456 VStep *= 1.1;
1457 V+=VStep;
1458 if(V*sv.VStep > sv.VStop*sv.VStep && V*sv.VStep < (sv.VStop + VStep - 1e-10*VStep)*sv.VStep)
1459 V=sv.VStop;
1460 }
1461 else // oh, diverged... reduce step and try again
1462 {
1463 if(++diverged_retry>=8)
1464 {
1465 gss_log.string_buf()<<"------> Too many failed steps, give up tring.\n\n\n";
1466 gss_log.record();
1467 break;
1468 }
1469 diverged_recovery();
1470 VStep /= 2.0;
1471 V-=VStep;
1472 gss_log.string_buf()<<"------> nonlinear solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1473 gss_log.record();
1474 continue;
1475 }
1476
1477 for(int j=0;j<sv.Electrode_Record.size();j++)
1478 {
1479 I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new();
1480 fprintf(fiv,"%lf\t%e\t",
1481 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential_new()/voltage_scale_V),
1482 double(I/current_scale_mA));
1483 }
1484 if(sv.Electrode_Record.size())
1485 {
1486 fprintf(fiv,"\n");
1487 fflush(fiv);
1488 }
1489
1490 }
1491 while(V*sv.VStep < (sv.VStop + 0.5*VStep)*sv.VStep);
1492 }
1493
1494 if(sv.Electrode_IScan!=-1)
1495 {
1496 bc.Set_electrode_type(sv.Electrode_IScan_Name.c_str(),CurrentBC);
1497 PetscScalar I = sv.IStart;
1498 PetscScalar IStep = sv.IStep;
1499 int diverged_retry=0;
1500 probe_open(DCSWEEP_ISCAN);
1501 do
1502 {
1503 its = 0;
1504 gss_log.string_buf()<<"DC Scan: I("<<sv.Electrode_IScan_Name<<") = "<<I/current_scale_mA<<" mA"<<"\n";
1505 gss_log.record();
1506 bc.Set_Iapp(sv.Electrode_IScan_Name.c_str(),I);
1507 solver_pre_compute();
1508 SNESSolve(snes,PETSC_NULL,x);
1509 SNESGetConvergedReason(snes,&reason);
1510
1511 if(reason>0 && norm==norm) //ok, converged.
1512 {
1513 diverged_retry=0;
1514 solution_update();
1515 probe(DCSWEEP_ISCAN,I);
1516 if(fabs(IStep) < fabs(sv.IStep))
1517 IStep *= 1.1;
1518 I+=IStep;
1519 if(I*sv.IStep > sv.IStop*sv.IStep && I*sv.IStep < (sv.IStop + IStep - 1e-10*IStep)*sv.IStep)
1520 I=sv.IStop;
1521 }
1522 else // oh, diverged... reduce step and try again
1523 {
1524 if(++diverged_retry>=8) //failed 8 times, stop tring
1525 {
1526 gss_log.string_buf()<<"------> Too many failed steps, give up tring.\n\n\n";
1527 gss_log.record();
1528 break;
1529 }
1530 diverged_recovery();
1531 IStep /= 2.0;
1532 I-=IStep;
1533 gss_log.string_buf()<<"------> nonlinear solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1534 gss_log.record();
1535 continue;
1536 }
1537
1538 for(int j=0;j<sv.Electrode_Record.size();j++)
1539 {
1540 fprintf(fiv,"%lf\t%e\t",
1541 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential_new()/voltage_scale_V),
1542 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new()/current_scale_mA));
1543 }
1544 if(sv.Electrode_Record.size())
1545 {
1546 fprintf(fiv,"\n");
1547 fflush(fiv);
1548 }
1549 }
1550 while(I*sv.IStep < (sv.IStop+0.5*IStep)*sv.IStep);
1551 }
1552
1553 if(!sv.IVFile.empty()) fclose(fiv);
1554 probe_close();
1555
1556 return 0;
1557
1558 }
1559
1560
1561
1562
1563 /* ----------------------------------------------------------------------------
1564 * QDDM_Solver_L1E::solve_transient: This function does transient simulation.
1565 */
1566 int QDDM_Solver_L1E::solve_transient(SolveDefine &sv)
1567 {
1568
1569 PetscScalar I=0;
1570 FILE *fiv;
1571
1572 PetscScalar current_scale_mA = scale_unit.s_mA;
1573 PetscScalar voltage_scale_V =scale_unit.s_volt;
1574 int diverged_retry=0;
1575
1576 clock = sv.TStart+sv.TStep;
1577 ODE_formula.TimeDependent = true;
1578 ODE_formula.BDF_Type = sv.BDF_Type;
1579 ODE_formula.dt = sv.TStep;
1580 if(sv.BDF_Type==BDF2)
1581 {
1582 ODE_formula.BDF2_restart = true;
1583 }
1584
1585 for(int i=0;i<zone_num;i++)
1586 {
1587 if(zonedata[i]->material_type == Semiconductor)
1588 {
1589 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[i]);
1590 pzonedata->HighFieldMobility= sv.HighFieldMobility;
1591 pzonedata->ImpactIonization = sv.ImpactIonization;
1592 pzonedata->IIType = sv.IIType;
1593 pzonedata->BandBandTunneling = sv.BandBandTunneling;
1594 pzonedata->IncompleteIonization = sv.IncompleteIonization;
1595 pzonedata->QuantumMechanical = sv.QuantumMechanical;
1596 pzonedata->Fermi = sv.Fermi;
1597 pzonedata->EJModel = sv.EJModel;
1598 pzonedata->QNFactor = sv.QNFactor;
1599 pzonedata->QPFactor = sv.QPFactor;
1600 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1601 for(int j=0;j<pzonedata->electrode.size();j++)
1602 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1603 {
1604 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1605 sv.Electrode_Record_Index.push_back(k);
1606 }
1607 }
1608 if(zonedata[i]->material_type == Insulator)
1609 {
1610 ISZone * pzonedata = dynamic_cast< ISZone * >(zonedata[i]);
1611 for(int k=0;k<sv.Electrode_Record_Name.size();k++)
1612 for(int j=0;j<pzonedata->electrode.size();j++)
1613 if(bc.is_electrode_label(pzonedata->electrode[j], sv.Electrode_Record_Name[k].c_str()))
1614 {
1615 sv.Electrode_Record.push_back(pzonedata->electrode[j]);
1616 sv.Electrode_Record_Index.push_back(k);
1617 }
1618 }
1619 }
1620
1621 gss_log.string_buf()<<"Transient compute from "<<sv.TStart<<" ps step "<<sv.TStep<<" ps to "<<sv.TStop<<" ps"<<"\n";
1622 gss_log.record();
1623
1624 //open record file
1625 if(!sv.IVFile.empty())
1626 {
1627 if(sv.IVFileAppend)
1628 {
1629 fiv=fopen(sv.IVFile.c_str(),"a");
1630 if(!fiv) fiv=stdout;
1631 }
1632 else
1633 {
1634 fiv=fopen(sv.IVFile.c_str(),"w");
1635 if(!fiv) fiv=stdout;
1636 fprintf(fiv,"#");
1637 fprintf(fiv,"time[ps] ");
1638 for(int j=0;j<sv.Electrode_Record.size();j++)
1639 fprintf(fiv,"Vp(%s)[V] I(%s)[mA] ",
1640 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str(),
1641 sv.Electrode_Record_Name[sv.Electrode_Record_Index[j]].c_str()
1642 );
1643 fprintf(fiv,"\n");
1644 }
1645 }
1646 else
1647 fiv = stdout;
1648
1649 probe_open(TRANSIENT);
1650
1651 // the main loop of transient solver.
1652 do
1653 {
1654 gss_log.string_buf()<<"t = "<<clock<<" ps"<<"\n";
1655 gss_log.record();
1656 ODE_formula.clock = clock;
1657 bc.Update_Vapp(clock);
1658 bc.Update_Iapp(clock);
1659 its = 0;
1660 optgen_update(clock);
1661 solver_pre_compute();
1662 SNESSolve(snes,PETSC_NULL,x);
1663 SNESGetConvergedReason(snes,&reason);
1664
1665 if(reason<0 || !(norm==norm))
1666 {
1667 if(++diverged_retry>=8) //failed 8 times, stop tring
1668 {
1669 gss_log.string_buf()<<"------> Too many failed steps, give up tring.\n\n\n";
1670 gss_log.record();
1671 break;
1672 }
1673 diverged_recovery();
1674 gss_log.string_buf()<<"------> nonlinear solver "<<SNESConvergedReasons[reason]<<", do recovery...\n\n\n";
1675 gss_log.record();
1676 ODE_formula.dt/=2.0;
1677 clock-=ODE_formula.dt;// reduce time step
1678 if(clock<sv.TStart) clock=sv.TStart;
1679 continue;
1680 }
1681 //ok, converged.
1682
1683 diverged_retry=0;
1684 solution_update();
1685 probe(TRANSIENT,clock);
1686 if(sv.Electrode_Record.size())
1687 {
1688 fprintf(fiv,"%e\t",double(clock));
1689 for(int j=0;j<sv.Electrode_Record.size();j++)
1690 {
1691 I = bc.Get_pointer(sv.Electrode_Record[j])->Get_Current_new();
1692 fprintf(fiv,"%lf\t%e\t",
1693 double(bc.Get_pointer(sv.Electrode_Record[j])->Get_Potential_new()/voltage_scale_V),
1694 double(I/current_scale_mA));
1695 }
1696 if(sv.Electrode_Record.size())
1697 {
1698 fprintf(fiv,"\n");
1699 fflush(fiv);
1700 }
1701 }
1702 if(fabs(ODE_formula.dt) < fabs(sv.TStep))
1703 ODE_formula.dt*=1.1;
1704 clock+=ODE_formula.dt;
1705 if(clock > sv.TStop && clock < (sv.TStop + ODE_formula.dt - 1e-10*ODE_formula.dt))
1706 {
1707 ODE_formula.dt -= clock - sv.TStop;
1708 clock=sv.TStop;
1709 }
1710
1711
1712 if(ODE_formula.BDF_Type==BDF2)
1713 ODE_formula.BDF2_restart = false;
1714
1715 }while(clock<sv.TStop+0.5*ODE_formula.dt);
1716
1717 //close record file
1718 if(!sv.IVFile.empty())
1719 fclose(fiv);
1720 probe_close();
1721 return 0;
1722
1723 }
1724
1725
1726
1727 /* ----------------------------------------------------------------------------
1728 * QDDM_Solver_L1E::solver_pre_compute: This function precompute some parameters before
1729 * each newton iteration.
1730 */
1731 void QDDM_Solver_L1E::solver_pre_compute()
1732 {
1733 for(int z=0;z<zone_num;z++)
1734 {
1735 if(zonedata[z]->material_type==Semiconductor)
1736 {
1737 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1738 for(int i=0;i<zone[z].davcell.size();i++)
1739 {
1740 const VoronoiCell* pcell = zone[z].davcell.GetPointer(i);
1741 SemiData *pfs = &(pzonedata->fs[i]);
1742 SemiAuxData *paux = &(pzonedata->aux[i]);
1743 pzonedata->mt->mapping(&zone[z].danode[i],paux,0);
1744 pzonedata->aux[i].mun = pzonedata->mt->mob->ElecMob(pfs->p,pfs->n,pfs->T,0.0,0.0,pfs->T);
1745 pzonedata->aux[i].mup = pzonedata->mt->mob->HoleMob(pfs->p,pfs->n,pfs->T,0.0,0.0,pfs->T);
1746 }
1747 }
1748 }
1749 }
1750
1751
1752 /* ----------------------------------------------------------------------------
1753 * QDDM_Solver_L1E::solver_post_compute: This function do post process after each
1754 * DC sweep or transient simulation
1755 */
1756 void QDDM_Solver_L1E::solver_post_compute()
1757 {
1758
1759 }
1760
1761
1762 /* ----------------------------------------------------------------------------
1763 * QDDM_Solver_L1E::solution_update: This function restore solution data from SNES
1764 */
1765 void QDDM_Solver_L1E::solution_update()
1766 {
1767 //------------------------------------------------------------
1768 PetscScalar *xx;
1769 VecGetArray(x,&xx);
1770 int offset=0;
1771 for(int z=0;z<zone_num;z++)
1772 {
1773 if(zonedata[z]->material_type==Semiconductor)
1774 {
1775 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1776 PetscScalar kb = pzonedata->mt->kb;
1777 PetscScalar e = pzonedata->mt->e;
1778 for(int i=0;i<zone[z].davcell.size();i++)
1779 {
1780 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
1781 pzonedata->fs[i].n_last = pzonedata->fs[i].n;
1782 pzonedata->fs[i].p_last = pzonedata->fs[i].p;
1783 pzonedata->fs[i].Eqc_last = pzonedata->fs[i].Eqc;
1784 pzonedata->fs[i].Eqv_last = pzonedata->fs[i].Eqv;
1785 pzonedata->fs[i].P = xx[offset+0];
1786 pzonedata->fs[i].n = xx[offset+1];
1787 pzonedata->fs[i].p = xx[offset+2];
1788 pzonedata->fs[i].Eqc = xx[offset+3];
1789 pzonedata->fs[i].Eqv = xx[offset+4];
1790
1791 pzonedata->mt->mapping(&pzonedata->pzone->danode[i],&pzonedata->aux[i],0);
1792 PetscScalar nie = pzonedata->mt->band->nie(pzonedata->fs[i].T);
1793 pzonedata->aux[i].phi_intrinsic = pzonedata->fs[i].P + pzonedata->aux[i].affinity +
1794 kb*pzonedata->fs[i].T/e*log(pzonedata->aux[i].Nc/nie);
1795 pzonedata->aux[i].phin = pzonedata->aux[i].phi_intrinsic - log(fabs(pzonedata->fs[i].n)/nie)*kb*pzonedata->fs[i].T/e;
1796 pzonedata->aux[i].phip = pzonedata->aux[i].phi_intrinsic + log(fabs(pzonedata->fs[i].p)/nie)*kb*pzonedata->fs[i].T/e;
1797
1798 offset += 5;
1799 }
1800 pzonedata->F1Q_efield_update(xx,zofs,bc,zonedata);
1801 for(int i=0;i<pzonedata->electrode.size();i++)
1802 {
1803 PetscScalar P = xx[offset];
1804 PetscScalar Pn = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
1805 PetscScalar C = bc.Get_pointer(pzonedata->electrode[i])->Get_C();
1806 PetscScalar I = bc.Get_pointer(pzonedata->electrode[i])->Get_Current_new();
1807 PetscScalar In= bc.Get_pointer(pzonedata->electrode[i])->Get_Current();
1808 bc.Get_pointer(pzonedata->electrode[i])->Set_Cap_Current(C*(P-Pn)/ODE_formula.dt);
1809 bc.Get_pointer(pzonedata->electrode[i])->Set_Current(I);
1810 bc.Get_pointer(pzonedata->electrode[i])->Set_Current_old(In);
1811 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential(P);
1812 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential_old(Pn);
1813 offset += 1;
1814 }
1815 }
1816 if(zonedata[z]->material_type==Insulator)
1817 {
1818 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1819 for(int i=0;i<zone[z].davcell.size();i++)
1820 {
1821 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
1822 pzonedata->fs[i].P = xx[offset];
1823 offset += 1;
1824 }
1825 pzonedata->F1_efield_update(xx,zofs,bc,zonedata);
1826 for(int i=0;i<pzonedata->electrode.size();i++)
1827 {
1828 PetscScalar P = xx[offset];
1829 PetscScalar Pn = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
1830 PetscScalar C = bc.Get_pointer(pzonedata->electrode[i])->Get_C();
1831 PetscScalar I = bc.Get_pointer(pzonedata->electrode[i])->Get_Current_new();
1832 PetscScalar In= bc.Get_pointer(pzonedata->electrode[i])->Get_Current();
1833 bc.Get_pointer(pzonedata->electrode[i])->Set_Cap_Current(C*(P-Pn)/ODE_formula.dt);
1834 bc.Get_pointer(pzonedata->electrode[i])->Set_Current(I);
1835 bc.Get_pointer(pzonedata->electrode[i])->Set_Current_old(In);
1836 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential(P);
1837 bc.Get_pointer(pzonedata->electrode[i])->Set_Potential_old(Pn);
1838 offset += 1;
1839 }
1840 }
1841 if(zonedata[z]->material_type==Conductor)
1842 {
1843 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1844 for(int i=0;i<zone[z].davcell.size();i++)
1845 {
1846 pzonedata->fs[i].P_last = pzonedata->fs[i].P;
1847 pzonedata->fs[i].P = xx[offset];
1848 offset += 1;
1849 }
1850 }
1851 }
1852 ODE_formula.dt_last = ODE_formula.dt;
1853
1854 VecRestoreArray(x,&xx);
1855 }
1856
1857 /* ----------------------------------------------------------------------------
1858 * QDDM_Solver_L1E::diverged_recovery: This function recovery latest solution data
1859 * if SNES diverged.
1860 */
1861 void QDDM_Solver_L1E::diverged_recovery()
1862 {
1863 //------------------------------------------------------------
1864 PetscScalar *xx;
1865 VecGetArray(x,&xx);
1866 int offset=0;
1867 for(int z=0;z<zone_num;z++)
1868 {
1869 if(zonedata[z]->material_type==Semiconductor)
1870 {
1871 SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
1872 for(int i=0;i<zone[z].davcell.size();i++)
1873 {
1874 xx[offset+0] = pzonedata->fs[i].P;
1875 xx[offset+1] = pzonedata->fs[i].n;
1876 xx[offset+2] = pzonedata->fs[i].p;
1877 xx[offset+3] = pzonedata->fs[i].Eqc;
1878 xx[offset+4] = pzonedata->fs[i].Eqv;
1879 offset += 5;
1880 }
1881 for(int i=0;i<pzonedata->electrode.size();i++)
1882 {
1883 PetscScalar P = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
1884 xx[offset] = P;
1885 offset += 1;
1886 }
1887 }
1888 if(zonedata[z]->material_type==Insulator)
1889 {
1890 ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
1891 for(int i=0;i<zone[z].davcell.size();i++)
1892 {
1893 xx[offset] = pzonedata->fs[i].P;
1894 offset += 1;
1895 }
1896 for(int i=0;i<pzonedata->electrode.size();i++)
1897 {
1898 PetscScalar P = bc.Get_pointer(pzonedata->electrode[i])->Get_Potential();
1899 xx[offset] = P;
1900 offset += 1;
1901 }
1902 }
1903 if(zonedata[z]->material_type==Conductor)
1904 {
1905 ElZone *pzonedata = dynamic_cast< ElZone * >(zonedata[z]);
1906 for(int i=0;i<zone[z].davcell.size();i++)
1907 {
1908 xx[offset] = pzonedata->fs[i].P;
1909 offset += 1;
1910 }
1911 }
1912 }
1913 VecRestoreArray(x,&xx);
1914 }
1915
1916
1917 /* ----------------------------------------------------------------------------
1918 * QDDM_Solver_L1E::destroy_solver: This function do destroy the nonlinear solver
1919 */
1920 int QDDM_Solver_L1E::destroy_solver(SolveDefine &sv)
1921 {
1922 // free work space
1923 N = 0;
1924 zofs.clear();
1925 VecDestroy(x);
1926 VecDestroy(r);
1927 MatDestroy(J);
1928 MatDestroy(JTmp);
1929 if(mf_flg)
1930 MatDestroy(JPrec);
1931 SNESDestroy(snes);
1932 return 0;
1933 }
General-purpose Semiconductor Simulator