bugfix Q_p calculation

This commit is contained in:
Markus Mueller 2020-05-11 15:48:03 +02:00
parent 67857cf983
commit 5bbe0f9cf7
1 changed files with 15 additions and 14 deletions

View File

@ -407,7 +407,7 @@ HICUMload(GENmodel *inModel, CKTcircuit *ckt)
double Qf_Vbiei, Qf_Vbici, Qf_Vciei, Qf_dT, Qf_dQ_pT, Qf_dick, Qf_dT_f0;
double Qr_Vbiei, Qr_Vbici, Qr_Vciei, Qr_dT, Qr_dQ_pT, Qr_dick, Qr_dT_f0;
double it_ditf, it_ditr;
duals::duald result_itf, result_itr, result_Qp, result_Qf, result_Qr, result_Q_bf, result_a_h; //intermediate variables when calling void dual functions
duals::duald result_itf, result_itr, result_Qp, result_Qf, result_Qr, result_Q_bf, result_a_h, result_Q_p; //intermediate variables when calling void dual functions
double Orci0_t,T_f0, Q_p, a_h;
double volatile Q_bf, Q_bf_Vbiei, Q_bf_Vbici, Q_bf_Vciei, Q_bf_dT, Q_bf_dick, Q_bf_dT_f0, Q_bf_dQ_pT;
double volatile Q_pT, Q_pT_dVbiei, Q_pT_dVbici, Q_pT_dT, Q_pT_dick, Q_pT_dT_f0, Q_pT_dQ_0, Q_pT_dVciei;
@ -547,7 +547,7 @@ HICUMload(GENmodel *inModel, CKTcircuit *ckt)
// Q_fC, Q_CT: actual and ICCR (weighted) hole charge
// T_fC, T_cT: actual and ICCR (weighted) transit time
// Derivative dfCT_ditf not properly implemented yet
// feenableexcept(FE_INVALID | FE_OVERFLOW); //debuger catches NANS
feenableexcept(FE_INVALID | FE_OVERFLOW); //debuger catches NANS
std::function<void (duals::duald, duals::duald, duals::duald, double, duals::duald*, duals::duald*, duals::duald*, duals::duald*)> HICQFC = [&](duals::duald T, duals::duald Ix, duals::duald I_CK, double FFT_pcS, duals::duald * Q_fC, duals::duald * Q_CT, duals::duald * T_fC, duals::duald * T_cT)
{
@ -896,24 +896,24 @@ HICUMload(GENmodel *inModel, CKTcircuit *ckt)
*Qr = Tr*(*itr);
};
std::function<void (duals::duald, duals::duald, duals::duald, duals::duald, duals::duald, duals::duald, duals::duald*, duals::duald*, duals::duald*, duals::duald*, duals::duald*, duals::duald*)> calc_it_initial = [&](duals::duald T, duals::duald Vbiei, duals::duald Vbici, duals::duald Q_0, duals::duald T_f0, duals::duald ick, duals::duald *itf, duals::duald *itr, duals::duald *Qf, duals::duald *Qr, duals::duald *Q_bf, duals::duald *a_h){
std::function<void (duals::duald, duals::duald, duals::duald, duals::duald, duals::duald, duals::duald, duals::duald*, duals::duald*, duals::duald*, duals::duald*, duals::duald*, duals::duald*, duals::duald*)> calc_it_initial = [&](duals::duald T, duals::duald Vbiei, duals::duald Vbici, duals::duald Q_0, duals::duald T_f0, duals::duald ick, duals::duald *itf, duals::duald *itr, duals::duald *Qf, duals::duald *Qr, duals::duald *Q_bf, duals::duald *a_h, duals::duald *Q_p){
// given T,Q_pT, ick, T_f0, Tr, Vbiei, Vbici -> calculate itf, itr, Qf, Qr
duals::duald VT, VT_f, i_0f, i_0r, I_Tf1, Tf, Q_fT, T_fT, A, Q_p;
duals::duald VT, VT_f, i_0f, i_0r, I_Tf1, Tf, Q_fT, T_fT, A;
VT = CONSTboltz * T / CHARGE;
VT_f = model->HICUMmcf*VT;
i_0f = here->HICUMc10_t * exp(Vbiei/VT_f);
i_0r = here->HICUMc10_t * exp(Vbici/VT);
Q_p = Q_0;
*Q_p = Q_0;
if (T_f0 > 0.0 || Tr > 0.0) {
A = 0.5*Q_0;
Q_p = A+sqrt(A*A+T_f0*i_0f+Tr*i_0r);
*Q_p = A+sqrt(A*A+T_f0*i_0f+Tr*i_0r);
}
I_Tf1 =i_0f/Q_p;
I_Tf1 =i_0f/(*Q_p);
*a_h = Oich*I_Tf1;
*itf = I_Tf1*(1.0+*a_h);
*itr = i_0r/Q_p;
*itr = i_0r/(*Q_p);
//Initial formulation of forward transit time, diffusion, GICCR and excess b-c charge
*Q_bf = 0.0;
@ -1633,13 +1633,14 @@ HICUMload(GENmodel *inModel, CKTcircuit *ckt)
//Initial formulation of forward and reverse component of transfer current
//begin initial transfer current calculations -> itf, itr, Qf, Qr------------
calc_it_initial(here->HICUMtemp+1_e, Vbiei , Vbici , Q_0 , T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h);
calc_it_initial(here->HICUMtemp+1_e, Vbiei , Vbici , Q_0 , T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h, &result_Q_p);
itf = result_itf.rpart();
itr = result_itr.rpart();
Qf = result_Qf.rpart();
Qr = result_Qr.rpart();
Q_bf = result_Q_bf.rpart();
a_h = result_a_h.rpart(); //needed to check if newton iteration needed
Q_p = result_Q_p.rpart(); //needed to check if newton iteration needed
itf_dT = result_itf.dpart();
itr_dT = result_itr.dpart();
Qf_dT = result_Qf.dpart();
@ -1647,35 +1648,35 @@ HICUMload(GENmodel *inModel, CKTcircuit *ckt)
Q_bf_dT= result_Q_bf.dpart();
if (!(Qf > RTOLC*Q_p || a_h > RTOLC)) { // in this case the newon is not run and the derivatives of the initial solution are needed
calc_it_initial(here->HICUMtemp , Vbiei+1_e, Vbici , Q_0 , T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h);
calc_it_initial(here->HICUMtemp , Vbiei+1_e, Vbici , Q_0 , T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h, &result_Q_p);
itf_Vbiei = result_itf.dpart();
itr_Vbiei = result_itr.dpart();
Qf_Vbiei = result_Qf.dpart();
Qr_Vbiei = result_Qr.dpart();
Q_bf_Vbiei = result_Q_bf.dpart();
calc_it_initial(here->HICUMtemp , Vbiei , Vbici+1_e, Q_0 , T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h);
calc_it_initial(here->HICUMtemp , Vbiei , Vbici+1_e, Q_0 , T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h, &result_Q_p);
itf_Vbici = result_itf.dpart();
itr_Vbici = result_itr.dpart();
Qf_Vbici = result_Qf.dpart();
Qr_Vbici = result_Qr.dpart();
Q_bf_Vbici= result_Q_bf.dpart();
calc_it_initial(here->HICUMtemp , Vbiei , Vbici , Q_0+1_e, T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h);
calc_it_initial(here->HICUMtemp , Vbiei , Vbici , Q_0+1_e, T_f0 , ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h, &result_Q_p);
itf_dQ_pT = result_itf.dpart(); //actual Q_0, but we can spare us this new variables
itr_dQ_pT = result_itr.dpart();
Qf_dQ_pT = result_Qf.dpart();
Qr_dQ_pT = result_Qr.dpart();
Q_bf_dQ_pT= result_Q_bf.dpart();
calc_it_initial(here->HICUMtemp , Vbiei , Vbici , Q_0 , T_f0+1_e, ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h);
calc_it_initial(here->HICUMtemp , Vbiei , Vbici , Q_0 , T_f0+1_e, ick , &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h, &result_Q_p);
itf_dT_f0 = result_itf.dpart();
itr_dT_f0 = result_itr.dpart();
Qf_dT_f0 = result_Qf.dpart();
Qr_dT_f0 = result_Qr.dpart();
Q_bf_dT_f0= result_Q_bf.dpart();
calc_it_initial(here->HICUMtemp , Vbiei , Vbici , Q_0 , T_f0 , ick+1_e, &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h);
calc_it_initial(here->HICUMtemp , Vbiei , Vbici , Q_0 , T_f0 , ick+1_e, &result_itf, &result_itr, &result_Qf, &result_Qr, &result_Q_bf, &result_a_h, &result_Q_p);
itf_dick = result_itf.dpart();
itr_dick = result_itr.dpart();
Qf_dick = result_Qf.dpart();