initial checkin, based on GSS 0.46 CVS
[gss-tcad.git] / src / solver / mix2 / insulequ2mix.cc
blobeb36adecfcb213e60c0992935dec97e4999ab53a
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, 2005 */
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 /*****************************************************************************/
21 #include "zonedata.h"
23 //-----------------------------------------------------------------------------
24 // special process of electrode bcs for mix type simulation
25 //-----------------------------------------------------------------------------
27 void ISZone::F2_mix_ddm_gatebc_segment(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs,DABC &bc )
29 const VoronoiCell* pcell = &pzone->davcell[i];
30 GateBC *pbc = dynamic_cast<GateBC * >(bc.Get_pointer(pcell->bc_index-1));
31 PetscScalar Ti = x[zofs[zone_index]+2*i+1];
32 PetscScalar grad_T=0;
33 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
34 for(int j=0;j<pcell->nb_num;j++)
36 int nb = pcell->nb_array[j];
37 PetscScalar Tj = x[zofs[zone_index]+2*nb+1]; //potential of nb node
38 PetscScalar T_mid = 0.5*(Tj+Ti);
39 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
40 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
41 grad_T += kapa*pcell->elen[j]*dTdx_mid/pcell->area;
42 if(j==0||j==pcell->nb_num-1)
44 PetscScalar h = pbc->heat_transfer;
45 PetscScalar r = h*pbc->T_external;
46 grad_T += 0.5*pcell->ilen[j]*(r-0.25*h*(3*Ti+Tj))/pcell->area;
50 f[zofs[zone_index]+2*i+0] = x[zofs[zone_index]+2*i+0] - pbc->Vapp+ pbc->WorkFunction;
51 f[zofs[zone_index]+2*i+1] = grad_T;
53 if(ODE_F.TimeDependent==true)
55 //second order
56 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
58 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
59 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
60 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
61 f[zofs[zone_index]+2*i+1] += -aux[i].density*HeatCapacity*Tt/(ODE_F.dt_last+ODE_F.dt);
63 else //first order
65 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
66 f[zofs[zone_index]+2*i+1] += -aux[i].density*HeatCapacity*(Ti-fs[i].T)/ODE_F.dt;
70 //just for fill the extra equ entry
71 int equ_num = 2;
72 int size = pzone->davcell.size();
73 int gate_equ;
74 for(int j=0;j<electrode.size();j++)
76 if(electrode[j]==pcell->bc_index-1) {gate_equ=j;break;}
78 f[zofs[zone_index]+equ_num*size+gate_equ]=0.0;
82 void ISZone::F2_mix_ddm_gatebc_interface(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs,DABC &bc,
83 ElZone *pz, int n)
85 const VoronoiCell* pcell = &pzone->davcell[i];
86 GateBC *pbc = dynamic_cast<GateBC * >(bc.Get_pointer(pcell->bc_index-1));
87 PetscScalar Ti = x[zofs[zone_index]+2*i+1];
88 PetscScalar grad_T=0;
89 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
90 for(int j=0;j<pcell->nb_num;j++)
92 int nb = pcell->nb_array[j];
93 PetscScalar Tj = x[zofs[zone_index]+2*nb+1]; //potential of nb node
94 PetscScalar T_mid = 0.5*(Tj+Ti);
95 PetscScalar dTdx_mid = (Tj-Ti)/pcell->ilen[j];
96 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
97 grad_T += kapa*pcell->elen[j]*dTdx_mid;
99 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
100 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
101 for(int j=0;j<ncell->nb_num;j++)
103 int nb = ncell->nb_array[j];
104 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1]; //potential of nb node
105 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
106 grad_T += kapa*ncell->elen[j]/ncell->ilen[j]*(Tj_n-Ti);
109 int equ_num = 2;
110 int size = pzone->davcell.size();
111 int gate_equ;
112 for(int j=0;j<electrode.size();j++)
114 if(electrode[j]==pcell->bc_index-1) {gate_equ=j;break;}
116 f[zofs[zone_index]+2*i+0] = x[zofs[zone_index]+2*i+0] - pbc->Vapp + pbc->WorkFunction;
117 f[zofs[zone_index]+2*i+1] = grad_T;
119 if(ODE_F.TimeDependent==true)
121 //second order
122 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
124 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
125 PetscScalar Tt = (2-r)/(1-r)*Ti-1.0/(r*(1-r))*fs[i].T+(1-r)/r*fs[i].T_last;
126 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
127 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
128 f[zofs[zone_index]+2*i+1] += -aux[i].density*HeatCapacity1*Tt/(ODE_F.dt_last+ODE_F.dt)*pcell->area
129 -pz->aux[n].density*HeatCapacity2*Tt/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
131 else //first order
133 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
134 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
135 f[zofs[zone_index]+2*i+1] += -aux[i].density*HeatCapacity1*(Ti-fs[i].T)/ODE_F.dt*pcell->area
136 -pz->aux[n].density*HeatCapacity2*(Ti-fs[i].T)/ODE_F.dt*ncell->area;
140 //just for fill the extra equ entry
141 f[zofs[zone_index]+equ_num*size+gate_equ]=0.0;
145 void ISZone::F2_mix_gate_electrode_current(int i,PetscScalar *x,PetscScalar *f, ODE_Formula &ODE_F, vector<int> & zofs, DABC &bc, PetscScalar DeviceDepth)
147 int bc_index = electrode[i];
148 GateBC *pbc = dynamic_cast <GateBC * > (bc.Get_pointer(bc_index));
149 PetscScalar current=0;
150 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
152 int node=bc[bc_index].psegment->node_array[j];
153 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
154 PetscScalar Vi = x[zofs[zone_index]+2*node+0]; //potential of node i
156 for(int k=0;k<pcell->nb_num;k++)
158 int nb = pcell->nb_array[k];
159 PetscScalar Vj = x[zofs[zone_index]+2*nb+0]; //potential of nb node
160 //displacement current
161 current += DeviceDepth*pcell->elen[k]*aux[node].eps*((Vi-Vj)-(fs[node].P-fs[nb].P))/pcell->ilen[k]/ODE_F.dt;
164 pbc->Set_Current_new(current);
167 //-------------------------------------------------------------------------------------------------
169 void ISZone::J2_mix_ddm_gatebc_segment(int i,PetscScalar *x,Mat *jac, ODE_Formula &ODE_F, vector<int> & zofs,DABC &bc )
171 const VoronoiCell* pcell = &pzone->davcell[i];
172 GateBC *pbc = dynamic_cast<GateBC * >(bc.Get_pointer(pcell->bc_index-1));
173 PetscScalar Ti = x[zofs[zone_index]+2*i+1];
174 PetscScalar d_grad_T_dTi=0;
175 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
176 for(int j=0;j<pcell->nb_num;j++)
178 int nb = pcell->nb_array[j];
180 PetscScalar Tj = x[zofs[zone_index]+2*nb+1]; //potential of nb node
181 PetscScalar T_mid = 0.5*(Tj+Ti);
182 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
183 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j]/pcell->area;
184 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j]/pcell->area;
185 if(j==0||j==pcell->nb_num-1)
187 PetscScalar h = pbc->heat_transfer;
188 d_grad_T_dTi += -0.5*pcell->ilen[j]*0.25*h*3/pcell->area;
189 d_grad_T_dTj += -0.5*pcell->ilen[j]*0.25*h/pcell->area;
191 MatSetValue(*jac,zofs[zone_index]+2*i+1,zofs[zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
193 MatSetValue(*jac,zofs[zone_index]+2*i+0,zofs[zone_index]+2*i+0,1.0,INSERT_VALUES);
195 if(ODE_F.TimeDependent==true)
197 //second order
198 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
200 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
201 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
202 d_grad_T_dTi += -aux[i].density*HeatCapacity*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt);
204 else //first order
206 PetscScalar HeatCapacity = mt->thermal->HeatCapacity(Ti);
207 d_grad_T_dTi += -aux[i].density*HeatCapacity/ODE_F.dt;
210 MatSetValue(*jac,zofs[zone_index]+2*i+1,zofs[zone_index]+2*i+1,d_grad_T_dTi,INSERT_VALUES);
213 //just for fill the matrix entry
214 int equ_num = 2;
215 int size = pzone->davcell.size();
216 int gate_equ;
217 for(int j=0;j<electrode.size();j++)
219 if(electrode[j]==pcell->bc_index-1) {gate_equ=j;break;}
221 MatSetValue(*jac,zofs[zone_index]+equ_num*size+gate_equ,zofs[zone_index]+equ_num*size+gate_equ,1.0,INSERT_VALUES);
225 void ISZone::J2_mix_ddm_gatebc_interface(int i,PetscScalar *x,Mat *jac, ODE_Formula &ODE_F, vector<int> & zofs,DABC &bc,
226 ElZone *pz, int n)
228 const VoronoiCell* pcell = &pzone->davcell[i];
229 GateBC *pbc = dynamic_cast<GateBC * >(bc.Get_pointer(pcell->bc_index-1));
230 PetscScalar Ti = x[zofs[zone_index]+2*i+1];
231 PetscScalar d_grad_T_dTi=0;
232 mt->mapping(&pzone->danode[i],&aux[i],ODE_F.clock);
233 for(int j=0;j<pcell->nb_num;j++)
235 int nb = pcell->nb_array[j];
236 PetscScalar Tj = x[zofs[zone_index]+2*nb+1]; //potential of nb node
237 PetscScalar T_mid = 0.5*(Tj+Ti);
238 PetscScalar kapa = mt->thermal->HeatConduction(T_mid);
239 PetscScalar d_grad_T_dTj = kapa*pcell->elen[j]/pcell->ilen[j];
240 d_grad_T_dTi += -kapa*pcell->elen[j]/pcell->ilen[j];
241 MatSetValue(*jac,zofs[zone_index]+2*i+1,zofs[zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
243 const VoronoiCell* ncell = pz->pzone->davcell.GetPointer(n);
244 pz->mt->mapping(&pz->pzone->danode[n],&pz->aux[n],ODE_F.clock);
245 for(int j=0;j<ncell->nb_num;j++)
247 int nb = ncell->nb_array[j];
248 PetscScalar Tj_n = x[zofs[pz->pzone->zone_index]+2*nb+1];
249 PetscScalar kapa = pz->mt->thermal->HeatConduction(0.5*(Ti+Tj_n));
250 d_grad_T_dTi += -kapa*ncell->elen[j]/ncell->ilen[j];
251 PetscScalar d_grad_T_dTj = kapa*ncell->elen[j]/ncell->ilen[j];
252 MatSetValue(*jac,zofs[zone_index]+2*i+1,zofs[pz->pzone->zone_index]+2*nb+1,d_grad_T_dTj,INSERT_VALUES);
255 MatSetValue(*jac,zofs[zone_index]+2*i+0,zofs[zone_index]+2*i+0,1.0,INSERT_VALUES);
257 if(ODE_F.TimeDependent==true)
259 //second order
260 if(ODE_F.BDF_Type==BDF2&&ODE_F.BDF2_restart==false)
262 PetscScalar r = ODE_F.dt_last/(ODE_F.dt_last+ODE_F.dt);
263 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
264 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
265 d_grad_T_dTi += -aux[i].density*HeatCapacity1*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*pcell->area
266 -pz->aux[n].density*HeatCapacity2*(2-r)/(1-r)/(ODE_F.dt_last+ODE_F.dt)*ncell->area;
268 else //first order
270 PetscScalar HeatCapacity1 = mt->thermal->HeatCapacity(Ti);
271 PetscScalar HeatCapacity2 = pz->mt->thermal->HeatCapacity(Ti);
272 d_grad_T_dTi += -aux[i].density*HeatCapacity1/ODE_F.dt*pcell->area
273 -pz->aux[n].density*HeatCapacity2/ODE_F.dt*ncell->area;
276 MatSetValue(*jac,zofs[zone_index]+2*i+1,zofs[zone_index]+2*i+1,d_grad_T_dTi,INSERT_VALUES);
278 //just for fill the matrix entry
279 int equ_num = 2;
280 int size = pzone->davcell.size();
281 int gate_equ;
282 for(int j=0;j<electrode.size();j++)
283 if(electrode[j]==pcell->bc_index-1)
284 {gate_equ=j;break;}
285 MatSetValue(*jac,zofs[zone_index]+equ_num*size+gate_equ,zofs[zone_index]+equ_num*size+gate_equ,1.0,INSERT_VALUES);
289 void ISZone::F2_mix_gate_electrode_Load(int bc_index, double & current, PetscScalar & pdI_pdV, Vec & pdI_pdw, Vec & pdF_pdV,PetscScalar *x, Mat *jac, ODE_Formula &ODE_F, vector<int> &zofs, DABC & bc, PetscScalar DeviceDepth)
291 GateBC *pbc = dynamic_cast <GateBC * > (bc.Get_pointer(bc_index));
292 pdI_pdV = 0;
293 PetscScalar e = mt->e;
294 VecZeroEntries(pdI_pdw);
295 VecZeroEntries(pdF_pdV);
296 PetscScalar * apdI_pdw;
297 PetscScalar * apdF_pdV;
298 VecGetArray(pdI_pdw,&apdI_pdw);
299 VecGetArray(pdF_pdV,&apdF_pdV);
300 current = pbc->Get_Current_new();
302 for(int j=0;j<bc[bc_index].psegment->node_array.size();j++)
304 int node=bc[bc_index].psegment->node_array[j];
305 const VoronoiCell* pcell = pzone->davcell.GetPointer(node);
306 for(int k=0;k<pcell->nb_num;k++)
308 int nb = pcell->nb_array[k];
309 //for displacement current
310 PetscScalar dJ_dVi = aux[node].eps/pcell->ilen[k]/ODE_F.dt*pcell->elen[k];
311 PetscScalar dJ_dVr = -aux[node].eps/pcell->ilen[k]/ODE_F.dt*pcell->elen[k];
312 apdI_pdw[zofs[zone_index]+2*node] += dJ_dVi*DeviceDepth;
313 apdI_pdw[zofs[zone_index]+2*nb] += dJ_dVr*DeviceDepth;
315 apdF_pdV[zofs[zone_index]+2*node] = 1.0;
316 pdI_pdV = 0;
318 VecRestoreArray(pdI_pdw,&apdI_pdw);
319 VecRestoreArray(pdF_pdV,&apdF_pdV);