- Working Ebers Moll model. That's a significant step ahead.
- Simple 2x2|RHS SPICE stamps now supported by two terminal devices.
  This was implicitly contained in the design, but set_mat now shows how
  a two-terminal device fits into a SPICE approach.
- Introduced direct solvers for net groups with 1 or 2 nets.
- Introduced specialized solvers for N=1,2,3,4,5 by using templates.
- nice performance increase for BJTs
This commit is contained in:
Couriersud 2014-01-18 13:58:15 +00:00
parent 34f2318c02
commit 5876abcfaf
12 changed files with 853 additions and 111 deletions

View File

@ -16,6 +16,7 @@ NETLIST_START(7400_astable)
NETDEV_SOLVER(Solver)
NETDEV_PARAM(Solver.FREQ, 48000)
NETDEV_PARAM(Solver.ACCURACY, 1e-7)
// astable NAND Multivibrator
NETDEV_R(R1, 1000)

View File

@ -18,7 +18,7 @@ NETLIST_START(bjt)
/* NPN - example */
NETDEV_QJT_SW(Q, "BC237B")
NETDEV_QBJT_SW(Q, "BC237B")
NETDEV_R(RB, 1000)
NETDEV_R(RC, 1000)
@ -30,7 +30,7 @@ NETLIST_START(bjt)
/* PNP - example */
NETDEV_QPNP(Q1, "BC556B")
NETDEV_QBJT_SW(Q1, "BC556B")
NETDEV_R(RB1, 1000)
NETDEV_R(RC1, 1000)
@ -40,7 +40,7 @@ NETLIST_START(bjt)
NET_C(RB1.2, Q1.B)
NET_C(Q1.E, V3)
//NETDEV_LOG(logB, Q1.B)
//NETDEV_LOG(logC, Q1.C)
NETDEV_LOG(logB, Q.B)
NETDEV_LOG(logC, Q.C)
NETLIST_END()

View File

@ -10,10 +10,11 @@ NETLIST_START(bjt)
/* Standard stuff */
NETDEV_CLOCK(clk)
NETDEV_PARAM(clk.FREQ, 1000) // 1000 Hz
NETDEV_PARAM(clk.FREQ, 10000) // 1000 Hz
NETDEV_SOLVER(Solver)
NETDEV_PARAM(Solver.FREQ, 48000)
NETDEV_PARAM(Solver.ACCURACY, 1e-4)
NETDEV_PARAM(Solver.ACCURACY, 1e-6)
NETDEV_PARAM(Solver.RESCHED_LOOPS, 50)
NETDEV_ANALOG_INPUT(V5, 5)
NETDEV_ANALOG_INPUT(V3, 3.5)
@ -30,7 +31,13 @@ NETLIST_START(bjt)
NET_C(RB.2, Q.B)
NET_C(Q.E, GND)
//NETDEV_LOG(logB, Q.B)
//NETDEV_LOG(logC, Q.C)
// put some load on Q.C
NETDEV_R(RCE, 150000)
NET_C(RCE.1, Q.C)
NET_C(RCE.2, GND)
NETDEV_LOG(logB, Q.B)
NETDEV_LOG(logC, Q.C)
NETLIST_END()

View File

@ -13,7 +13,9 @@ NETLIST_START(msx)
NETDEV_PARAM(clk.FREQ, 1000) // 1000 Hz
NETDEV_SOLVER(Solver)
NETDEV_PARAM(Solver.FREQ, 48000)
NETDEV_PARAM(Solver.ACCURACY, 1e-4)
NETDEV_PARAM(Solver.ACCURACY, 1e-6)
NETDEV_PARAM(Solver.CONVERG, 0.3)
NETDEV_PARAM(Solver.RESCHED_LOOPS, 60)
NETDEV_R(RAY8910, 2345) // Max Voltage

View File

@ -49,10 +49,16 @@ NETLIB_RESET(Q)
{
}
NETLIB_UPDATE(Q)
{
// netlist().solver()->schedule1();
}
// ----------------------------------------------------------------------------------------
// nld_QBJT_switch
// ----------------------------------------------------------------------------------------
#if USE_OLD_S
NETLIB_START(QBJT_switch)
{
NETLIB_NAME(Q)::start();
@ -106,19 +112,80 @@ NETLIB_START(QBJT_switch)
}
NETLIB_UPDATE(Q)
NETLIB_UPDATE_PARAM(QBJT_switch)
{
netlist().solver()->schedule1();
}
#else
NETLIB_START(QBJT_switch)
{
NETLIB_NAME(Q)::start();
register_terminal("B", m_RB.m_P);
register_terminal("E", m_RB.m_N);
register_terminal("C", m_RC.m_P);
register_terminal("_E1", m_RC.m_N);
register_terminal("_B1", m_BC_dummy.m_P);
register_terminal("_C1", m_BC_dummy.m_N);
connect(m_RB.m_N, m_RC.m_N);
connect(m_RB.m_P, m_BC_dummy.m_P);
connect(m_RC.m_P, m_BC_dummy.m_N);
save(NAME(m_state_on));
m_RB.set(NETLIST_GMIN, 0.0, 0.0);
m_RC.set(NETLIST_GMIN, 0.0, 0.0);
m_BC_dummy.set(NETLIST_GMIN, 0.0, 0.0);
m_state_on = 0;
{
double IS = m_model.model_value("IS", 1e-15);
double BF = m_model.model_value("BF", 100);
double NF = m_model.model_value("NF", 1);
//double VJE = m_model.dValue("VJE", 0.75);
set_qtype((m_model.model_type() == "NPN") ? BJT_NPN : BJT_PNP);
double alpha = BF / (1.0 + BF);
diode d(IS, NF);
// Assume 5mA Collector current for switch operation
m_V = d.V(0.005 / alpha);
/* Base current is 0.005 / beta
* as a rough estimate, we just scale the conductance down */
m_gB = d.gI(0.005 / alpha);
if (m_gB < NETLIST_GMIN)
m_gB = NETLIST_GMIN;
m_gC = d.gI(0.005); // very rough estimate
//printf("%f %f \n", m_V, m_gB);
}
}
NETLIB_UPDATE_PARAM(QBJT_switch)
{
}
#endif
// ----------------------------------------------------------------------------------------
// nld_Q - Ebers Moll
// ----------------------------------------------------------------------------------------
#define USE_OLD_B (0)
#if USE_OLD_B
NETLIB_START(QBJT_EB)
{
NETLIB_NAME(Q)::start();
@ -154,6 +221,7 @@ NETLIB_START(QBJT_EB)
//double VJE = m_model.dValue("VJE", 0.75);
set_qtype((m_model.model_type() == "NPN") ? BJT_NPN : BJT_PNP);
//printf("type %s\n", m_model.model_type().cstr());
m_alpha_f = BF / (1.0 + BF);
m_alpha_r = BR / (1.0 + BR);
@ -163,9 +231,58 @@ NETLIB_START(QBJT_EB)
}
}
NETLIB_RESET(QBJT_EB)
NETLIB_UPDATE(QBJT_EB)
{
NETLIB_NAME(Q)::reset();
#if !(USE_ALTERNATE_SCHEDULING)
netlist().solver()->schedule1();
#else
if (!m_D_BE.m_P.net().isRailNet())
m_D_BE.m_P.net().solve(); // Basis
else if (!m_D_BE.m_N.net().isRailNet())
m_D_BE.m_N.net().solve(); // Emitter
else
m_D_BC.m_N.net().solve(); // Collector
#endif
}
#else
NETLIB_START(QBJT_EB)
{
NETLIB_NAME(Q)::start();
register_terminal("E", m_D_EB.m_P); // Cathode
register_terminal("B", m_D_EB.m_N); // Anode
register_terminal("C", m_D_CB.m_P); // Cathode
register_terminal("_B1", m_D_CB.m_N); // Anode
register_terminal("_E1", m_D_EC.m_P);
register_terminal("_C1", m_D_EC.m_N);
connect(m_D_EB.m_P, m_D_EC.m_P);
connect(m_D_EB.m_N, m_D_CB.m_N);
connect(m_D_CB.m_P, m_D_EC.m_N);
m_gD_BE.save("m_D_BE", *this);
m_gD_BC.save("m_D_BC", *this);
{
double IS = m_model.model_value("IS", 1e-15);
double BF = m_model.model_value("BF", 100);
double NF = m_model.model_value("NF", 1);
double BR = m_model.model_value("BR", 1);
double NR = m_model.model_value("NR", 1);
//double VJE = m_model.dValue("VJE", 0.75);
set_qtype((m_model.model_type() == "NPN") ? BJT_NPN : BJT_PNP);
//printf("type %s\n", m_model.model_type().cstr());
m_alpha_f = BF / (1.0 + BF);
m_alpha_r = BR / (1.0 + BR);
m_gD_BE.set_param(IS / m_alpha_f, NF);
m_gD_BC.set_param(IS / m_alpha_r, NR);
}
}
NETLIB_UPDATE(QBJT_EB)
@ -173,11 +290,21 @@ NETLIB_UPDATE(QBJT_EB)
#if !(USE_ALTERNATE_SCHEDULING)
netlist().solver()->schedule1();
#else
m_D_BE.m_P.net().solve();
m_D_BE.m_N.net().solve();
m_D_BC.m_P.net().solve();
if (!m_D_EB.m_P.net().isRailNet())
m_D_EB.m_P.net().solve(); // Basis
else if (!m_D_EB.m_N.net().isRailNet())
m_D_EB.m_N.net().solve(); // Emitter
else
m_D_CB.m_N.net().solve(); // Collector
#endif
}
#endif
NETLIB_RESET(QBJT_EB)
{
NETLIB_NAME(Q)::reset();
}
NETLIB_UPDATE_PARAM(QBJT_EB)

View File

@ -86,12 +86,16 @@ private:
* E
*/
#define USE_OLD_S (0)
#if USE_OLD_S
class NETLIB_NAME(QBJT_switch) : public NETLIB_NAME(QBJT)
{
public:
ATTR_COLD NETLIB_NAME(QBJT_switch)()
: NETLIB_NAME(QBJT)(BJT_SWITCH), m_gB(NETLIST_GMIN), m_gC(NETLIST_GMIN), m_V(0.0), m_state_on(0) { }
NETLIB_UPDATEI()
{
double vE = INPANALOG(m_EV);
@ -114,6 +118,7 @@ public:
m_RB.set(gb, v, 0.0);
m_RC.set(gc, 0.0, 0.0);
m_state_on = new_state;
// Don't update m_RB - could cause an infinite loop
m_RB.update_dev();
m_RC.update_dev();
}
@ -138,11 +143,79 @@ protected:
private:
};
#else
class NETLIB_NAME(QBJT_switch) : public NETLIB_NAME(QBJT)
{
public:
ATTR_COLD NETLIB_NAME(QBJT_switch)()
: NETLIB_NAME(QBJT)(BJT_SWITCH),
m_RB(netlist_object_t::ANALOG),
m_RC(netlist_object_t::ANALOG),
m_BC_dummy(netlist_object_t::ANALOG),
m_gB(NETLIST_GMIN), m_gC(NETLIST_GMIN), m_V(0.0), m_state_on(0) { }
NETLIB_UPDATE_TERMINALS()
{
double m = (is_qtype( BJT_NPN) ? 1 : -1);
int new_state = (m_RB.deltaV() * m > m_V ) ? 1 : 0;
if (m_state_on ^ new_state)
{
double gb = m_gB;
double gc = m_gC;
double v = m_V * m;
if (!new_state )
{
// not conducting
gb = NETLIST_GMIN;
v = 0;
gc = NETLIST_GMIN;
}
m_RB.set(gb, v, 0.0);
m_RC.set(gc, 0.0, 0.0);
m_state_on = new_state;
}
}
ATTR_HOT ATTR_ALIGN void virtual update()
{
if (!m_RB.m_P.net().isRailNet())
m_RB.m_P.net().solve(); // Basis
else if (!m_RB.m_N.net().isRailNet())
m_RB.m_N.net().solve(); // Emitter
else if (!m_RC.m_P.net().isRailNet())
m_RB.m_P.net().solve(); // Collector
}
nld_twoterm m_RB;
nld_twoterm m_RC;
// FIXME: the matrix solvers should be devices so we can properly
// schedule them. This is a workaround and blows netgroup size
nld_twoterm m_BC_dummy;
protected:
ATTR_COLD virtual void start();
ATTR_HOT void update_param();
double m_gB; // base conductance / switch on
double m_gC; // collector conductance / switch on
double m_V; // internal voltage source
UINT8 m_state_on;
private:
};
#endif
// ----------------------------------------------------------------------------------------
// nld_QBJT_EB
// ----------------------------------------------------------------------------------------
#define USE_OLD_B (0)
#if USE_OLD_B
class NETLIB_NAME(QBJT_EB) : public NETLIB_NAME(QBJT)
{
public:
@ -187,5 +260,57 @@ protected:
private:
};
#else
class NETLIB_NAME(QBJT_EB) : public NETLIB_NAME(QBJT)
{
public:
ATTR_COLD NETLIB_NAME(QBJT_EB)()
: NETLIB_NAME(QBJT)(BJT_EB),
m_D_CB(netlist_object_t::ANALOG),
m_D_EB(netlist_object_t::ANALOG),
m_D_EC(netlist_object_t::ANALOG)
{ }
NETLIB_UPDATE_TERMINALS()
{
m_gD_BE.update_diode(-m_D_EB.deltaV());
m_gD_BC.update_diode(-m_D_CB.deltaV());
double gee = m_gD_BE.G();
double gcc = m_gD_BC.G();
double gec = m_alpha_r * gcc;
double gce = m_alpha_f * gee;
double sIe = -m_gD_BE.I() + m_alpha_r * m_gD_BC.I();
double sIc = m_alpha_f * m_gD_BE.I() - m_gD_BC.I();
double Ie = sIe + gee * m_gD_BE.Vd() - gec * m_gD_BC.Vd();
double Ic = sIc - gce * m_gD_BE.Vd() + gcc * m_gD_BC.Vd();
//double Ie = sIe + gee * -m_D_EB.deltaV() - gec * -m_D_CB.deltaV();
//double Ic = sIc - gce * -m_D_EB.deltaV() + gcc * -m_D_CB.deltaV();
//printf("EB %f sIe %f sIc %f\n", m_D_BE.deltaV(), sIe, sIc);
m_D_EB.set_mat(gee, gec - gee, gce - gee, gee - gec, Ie, -Ie);
m_D_CB.set_mat(gcc, gce - gcc, gec - gcc, gcc - gce, Ic, -Ic);
m_D_EC.set_mat( 0, -gec, -gce, 0, 0, 0);
}
protected:
ATTR_COLD virtual void start();
ATTR_COLD virtual void reset();
ATTR_HOT void update_param();
ATTR_HOT ATTR_ALIGN void virtual update();
netlist_generic_diode m_gD_BC;
netlist_generic_diode m_gD_BE;
nld_twoterm m_D_CB; // gcc, gce - gcc, gec - gcc, gcc - gce | Ic
nld_twoterm m_D_EB; // gee, gec - gee, gce - gee, gee - gec | Ie
nld_twoterm m_D_EC; // 0, -gec, -gcc, 0 | 0
double m_alpha_f;
double m_alpha_r;
private:
};
#endif
#endif /* NLD_BJT_H_ */

View File

@ -20,9 +20,11 @@
ATTR_COLD void netlist_matrix_solver_t::setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &aowner)
{
/* make sure we loop at least once */
m_resched = true;
m_owner = &aowner;
NL_VERBOSE_OUT(("New solver setup\n"));
for (netlist_net_t * const * pn = nets.first(); pn != NULL; pn = nets.next(pn))
{
NL_VERBOSE_OUT(("setting up net\n"));
@ -33,6 +35,7 @@ ATTR_COLD void netlist_matrix_solver_t::setup(netlist_net_t::list_t &nets, NETLI
for (netlist_core_terminal_t *p = (*pn)->m_head; p != NULL; p = p->m_update_list_next)
{
NL_VERBOSE_OUT(("%s %s %d\n", p->name().cstr(), (*pn)->name().cstr(), (int) (*pn)->isRailNet()));
switch (p->type())
{
case netlist_terminal_t::TERMINAL:
@ -45,14 +48,21 @@ ATTR_COLD void netlist_matrix_solver_t::setup(netlist_net_t::list_t &nets, NETLI
case netlist_device_t::BJT_EB:
case netlist_device_t::DIODE:
//case netlist_device_t::VCVS:
//case netlist_device_t::BJT_SWITCH:
case netlist_device_t::BJT_SWITCH:
NL_VERBOSE_OUT(("found BJT/Diode\n"));
if (!m_dynamic.contains(&p->netdev()))
m_dynamic.add(&p->netdev());
break;
default:
break;
}
(*pn)->m_terms.add(static_cast<netlist_terminal_t *>(p));
{
netlist_terminal_t *pterm = static_cast<netlist_terminal_t *>(p);
if (pterm->m_otherterm->net().isRailNet())
(*pn)->m_rails.add(pterm);
else
(*pn)->m_terms.add(pterm);
}
NL_VERBOSE_OUT(("Added terminal\n"));
break;
case netlist_terminal_t::INPUT:
@ -65,17 +75,11 @@ ATTR_COLD void netlist_matrix_solver_t::setup(netlist_net_t::list_t &nets, NETLI
break;
}
}
NL_VERBOSE_OUT(("added net with %d populated connections (%d railnets)\n", (*pn)->m_terms.count(), (*pn)->m_rails.count()));
}
}
ATTR_HOT inline void netlist_matrix_solver_t::step(const netlist_time delta)
{
const double dd = delta.as_double();
for (netlist_core_device_t * const *p = m_steps.first(); p != NULL; p = m_steps.next(p))
(*p)->step_time(dd);
}
ATTR_HOT inline void netlist_matrix_solver_t::update_inputs()
ATTR_HOT ATTR_HOT inline void netlist_matrix_solver_t::update_inputs()
{
for (netlist_core_terminal_t * const *p = m_inps.first(); p != NULL; p = m_inps.next(p))
{
@ -91,7 +95,320 @@ ATTR_HOT inline void netlist_matrix_solver_t::update_inputs()
}
ATTR_HOT inline int netlist_matrix_solver_t::solve_non_dynamic()
ATTR_HOT ATTR_HOT inline void netlist_matrix_solver_t::update_dynamic()
{
/* update all non-linear devices */
for (netlist_core_device_t * const *p = m_dynamic.first(); p != NULL; p = m_dynamic.next(p))
switch ((*p)->family())
{
case netlist_device_t::DIODE:
static_cast<NETLIB_NAME(D) *>((*p))->update_terminals();
break;
default:
(*p)->update_terminals();
break;
}
}
void netlist_matrix_solver_t::schedule()
{
if (!solve())
{
// NL_VERBOSE_OUT(("update_inputs\n");
update_inputs();
}
else
{
//NL_VERBOSE_OUT(("resched\n");
if (m_owner != NULL)
this->m_owner->schedule1();
}
//solve();
// update_inputs();
}
ATTR_COLD void netlist_matrix_solver_t::reset()
{
m_last_step = netlist_time::zero;
}
ATTR_HOT inline void netlist_matrix_solver_t::step(const netlist_time delta)
{
const double dd = delta.as_double();
for (int k=0; k < m_steps.count(); k++)
m_steps[k]->step_time(dd);
}
// ----------------------------------------------------------------------------------------
// netlist_matrix_solver - Direct base
// ----------------------------------------------------------------------------------------
template <int m_N, int _storage_N>
ATTR_COLD int netlist_matrix_solver_directb_t<m_N, _storage_N>::get_net_idx(netlist_net_t *net)
{
for (int k = 0; k < N(); k++)
if (m_nets[k] == net)
return k;
return -1;
}
template <int m_N, int _storage_N>
ATTR_COLD void netlist_matrix_solver_directb_t<m_N, _storage_N>::setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner)
{
netlist_matrix_solver_t::setup(nets, owner);
m_term_num = 0;
m_rail_start = 0;
for (int k = 0; k < N(); k++)
{
netlist_net_t *net = m_nets[k];
const netlist_net_t::terminal_list_t &terms = net->m_terms;
for (int i = 0; i < terms.count(); i++)
{
m_terms[m_term_num].net_this = k;
int ot = get_net_idx(&terms[i]->m_otherterm->net());
m_terms[m_term_num].net_other = ot;
m_terms[m_term_num].term = terms[i];
if (ot>=0)
m_term_num++;
}
}
m_rail_start = m_term_num;
for (int k = 0; k < N(); k++)
{
netlist_net_t *net = m_nets[k];
const netlist_net_t::terminal_list_t &terms = net->m_terms;
const netlist_net_t::terminal_list_t &rails = net->m_rails;
for (int i = 0; i < terms.count(); i++)
{
m_terms[m_term_num].net_this = k;
int ot = get_net_idx(&terms[i]->m_otherterm->net());
m_terms[m_term_num].net_other = ot;
m_terms[m_term_num].term = terms[i];
if (ot<0)
m_term_num++;
}
for (int i = 0; i < rails.count(); i++)
{
m_terms[m_term_num].net_this = k;
m_terms[m_term_num].net_other = -1; //get_net_idx(&rails[i]->m_otherterm->net());
m_terms[m_term_num].term = rails[i];
m_term_num++;
}
}
}
template <int m_N, int _storage_N>
ATTR_HOT inline void netlist_matrix_solver_directb_t<m_N, _storage_N>::build_LE(double (* RESTRICT m_A)[_storage_N], double (* RESTRICT m_RHS))
{
#if 0
for (int i = 0; i < m_term_num; i++)
{
terms_t &t = m_terms[i];
m_RHS[t.net_this] += t.term->m_Idr;
m_A[t.net_this][t.net_this] += t.term->m_gt;
if (t.net_other >= 0)
{
//m_A[t.net_other][t.net_other] += t.term->m_otherterm->m_gt;
m_A[t.net_this][t.net_other] += -t.term->m_go;
//m_A[t.net_other][t.net_this] += -t.term->m_otherterm->m_go;
}
else
m_RHS[t.net_this] += t.term->m_go * t.term->m_otherterm->net().Q_Analog();
}
#else
for (int i = 0; i < m_rail_start; i++)
{
terms_t &t = m_terms[i];
m_RHS[t.net_this] += t.term->m_Idr;
m_A[t.net_this][t.net_this] += t.term->m_gt;
m_A[t.net_this][t.net_other] += -t.term->m_go;
}
for (int i = m_rail_start; i < m_term_num; i++)
{
terms_t &t = m_terms[i];
m_RHS[t.net_this] += t.term->m_Idr;
m_A[t.net_this][t.net_this] += t.term->m_gt;
m_RHS[t.net_this] += t.term->m_go * t.term->m_otherterm->net().Q_Analog();
}
#endif
}
// ----------------------------------------------------------------------------------------
// netlist_matrix_solver - Direct1
// ----------------------------------------------------------------------------------------
ATTR_HOT inline int netlist_matrix_solver_direct1_t::solve_non_dynamic()
{
#if 1
ATTR_UNUSED netlist_net_t *last_resched_net = NULL;
/* over-relaxation not really works on these matrices */
//const double w = 1.0; //2.0 / (1.0 + sin(3.14159 / (m_nets.count()+1)));
//const double w1 = 1.0 - w;
//NL_VERBOSE_OUT(("%d %d\n", N(), 1);
double gtot_t = 0.0;
double RHS_t = 0.0;
netlist_net_t *net = m_nets[0];
const netlist_net_t::terminal_list_t &rails = net->m_rails;
int rail_count = rails.count();
for (int i = 0; i < rail_count; i++)
{
gtot_t += rails[i]->m_gt;
RHS_t += rails[i]->m_Idr;
RHS_t += rails[i]->m_go * rails[i]->m_otherterm->net().Q_Analog();
}
double iIdr = RHS_t;
double new_val = iIdr / gtot_t;
#else
netlist_net_t *net = m_nets[0];
double m_A[1][1] = { {0.0} };
double m_RHS[1] = { 0.0 };
build_LE(m_A, m_RHS);
//NL_VERBOSE_OUT(("%f %f\n", new_val, m_RHS[0] / m_A[0][0]);
double new_val = m_RHS[0] / m_A[0][0];
#endif
double e = (new_val - net->m_cur.Analog);
double cerr = e * e;
net->m_cur.Analog = net->m_new.Analog = new_val;
if (is_dynamic() && (cerr > m_accuracy * m_accuracy))
{
return 2;
}
else
return 1;
}
ATTR_HOT inline bool netlist_matrix_solver_direct1_t::solve()
{
int resched_cnt = 0;
ATTR_UNUSED netlist_net_t *last_resched_net = NULL;
if (USE_ALTERNATE_SCHEDULING)
{
netlist_time now = owner().netlist().time();
netlist_time delta = now - m_last_step;
if (delta >= netlist_time::from_nsec(5)) // always update capacitors
{
NL_VERBOSE_OUT(("Step!\n"));
/* update all terminals for new time step */
m_last_step = now;
step(delta);
}
}
if (is_dynamic())
{
int this_resched;
do
{
update_dynamic();
this_resched = solve_non_dynamic();
resched_cnt += this_resched;
} while (this_resched > 1 && resched_cnt < m_resched_loops);
}
else
{
resched_cnt = solve_non_dynamic();
}
return (resched_cnt >= m_resched_loops);
}
// ----------------------------------------------------------------------------------------
// netlist_matrix_solver - Direct2
// ----------------------------------------------------------------------------------------
ATTR_HOT int netlist_matrix_solver_direct2_t::solve_non_dynamic()
{
double m_A[2][2] = { { 0.0 } };
double m_RHS[2] = { 0.0 };
build_LE(m_A, m_RHS);
//NL_VERBOSE_OUT(("%f %f\n", new_val, m_RHS[0] / m_A[0][0]);
const double a = m_A[0][0];
const double b = m_A[0][1];
const double c = m_A[1][0];
const double d = m_A[1][1];
const double new_val1 = a / (a*d - b*c) * (m_RHS[1] - c / a * m_RHS[0]);
const double new_val0 = (m_RHS[0] - b * new_val1) / a;
double e = (new_val0 - m_nets[0]->m_cur.Analog);
double cerr = e * e;
m_nets[0]->m_cur.Analog = m_nets[0]->m_new.Analog = new_val0;
e = (new_val1 - m_nets[1]->m_cur.Analog);
cerr += e * e;
m_nets[1]->m_cur.Analog = m_nets[1]->m_new.Analog = new_val1;
if (is_dynamic() && (cerr / N() > m_accuracy * m_accuracy))
{
return 2;
}
else
return 1;
}
ATTR_HOT inline bool netlist_matrix_solver_direct2_t::solve()
{
int resched_cnt = 0;
ATTR_UNUSED netlist_net_t *last_resched_net = NULL;
if (USE_ALTERNATE_SCHEDULING)
{
netlist_time now = owner().netlist().time();
netlist_time delta = now - m_last_step;
if (delta >= netlist_time::from_nsec(5)) // always update capacitors
{
NL_VERBOSE_OUT(("Step!\n"));
/* update all terminals for new time step */
m_last_step = now;
step(delta);
}
}
if (is_dynamic())
{
int this_resched;
do
{
update_dynamic();
this_resched = solve_non_dynamic();
resched_cnt += this_resched;
} while (this_resched > 1 && resched_cnt < m_resched_loops);
}
else
{
resched_cnt = solve_non_dynamic();
}
return (resched_cnt >= m_resched_loops);
}
// ----------------------------------------------------------------------------------------
// netlist_matrix_solver - Gauss - Seidel
// ----------------------------------------------------------------------------------------
template <int m_N, int _storage_N>
ATTR_HOT inline int netlist_matrix_solver_gauss_seidel_t<m_N, _storage_N>::solve_non_dynamic()
{
bool resched = false;
@ -99,45 +416,86 @@ ATTR_HOT inline int netlist_matrix_solver_t::solve_non_dynamic()
ATTR_UNUSED netlist_net_t *last_resched_net = NULL;
/* over-relaxation not really works on these matrices */
const double w = 1.0; //2.0 / (1.0 + sin(3.14159 / (m_nets.count()+1)));
const double w1 = 1.0 - w;
//const double w = 1.0; //2.0 / (1.0 + sin(3.14159 / (m_nets.count()+1)));
//const double w1 = 1.0 - w;
//printf("%f %d\n", w, m_nets.count());
double w[_storage_N];
double one_m_w[_storage_N];
double RHS[_storage_N];
for (int k = 0; k < N(); k++)
{
double gtot_t = 0.0;
double gabs_t = 0.0;
double RHS_t = 0.0;
netlist_net_t *net = m_nets[k];
const netlist_net_t::terminal_list_t &terms = net->m_terms;
const netlist_net_t::terminal_list_t &rails = net->m_rails;
const int term_count = terms.count();
const int rail_count = rails.count();
for (int i = 0; i < rail_count; i++)
{
gtot_t += rails[i]->m_gt;
gabs_t += fabs(rails[i]->m_go);
RHS_t += rails[i]->m_Idr;
RHS_t += rails[i]->m_go * rails[i]->m_otherterm->net().Q_Analog();
}
for (int i = 0; i < term_count; i++)
{
gtot_t += terms[i]->m_gt;
gabs_t += fabs(terms[i]->m_go);
RHS_t += terms[i]->m_Idr;
}
gabs_t *= m_convergence_factor;
if (gabs_t > gtot_t)
{
// Actually 1.0 / g_tot * g_tot / (gtot_t + gabs_t)
w[k] = 1.0 / (gtot_t + gabs_t);
one_m_w[k] = gabs_t / (gtot_t + gabs_t);
}
else
{
w[k] = 1.0 / gtot_t;
one_m_w[k] = 0.0;
}
RHS[k] = RHS_t;
}
//NL_VERBOSE_OUT(("%f %d\n", w, m_nets.count());
do {
resched = false;
double cerr = 0.0;
for (netlist_net_t * const *pn = m_nets.first(); pn != NULL; pn = m_nets.next(pn))
for (int k = 0; k < N(); k++)
{
netlist_net_t *net = *pn;
netlist_net_t *net = m_nets[k];
const netlist_net_t::terminal_list_t &terms = net->m_terms;
const int term_count = terms.count();
double gtot = 0;
double gabs = 0;
double iIdr = 0;
double new_val;
double iIdr = RHS[k];
for (int i = 0; i < terms.count(); i++)
for (int i = 0; i < term_count; i++)
{
gtot += terms[i]->m_gt;
gabs += fabs(terms[i]->m_go);
iIdr += terms[i]->m_Idr + terms[i]->m_go * terms[i]->m_otherterm->net().Q_Analog();
iIdr += terms[i]->m_go * terms[i]->m_otherterm->net().Q_Analog();
}
gabs *= m_convergence_factor;
if (gabs > gtot)
new_val = (net->m_cur.Analog * gabs + iIdr) / (gtot + gabs);
else
new_val = w1 * net->m_cur.Analog + w * iIdr / gtot;
if (fabs(new_val - net->m_cur.Analog) > m_accuracy)
{
resched = true;
last_resched_net = net;
}
//double new_val = (net->m_cur.Analog * gabs[k] + iIdr) / (gtot[k]);
double new_val = net->m_cur.Analog * one_m_w[k] + iIdr * w[k];
double e = (new_val - net->m_cur.Analog);
cerr += e * e;
net->m_cur.Analog = net->m_new.Analog = new_val;
NL_VERBOSE_OUT(("Info: %d\n", pn->object()->m_num_cons));
//NL_VERBOSE_OUT(("New: %lld %f %f\n", netlist().time().as_raw(), netlist().time().as_double(), new_val));
}
if (resched || cerr / m_nets.count() > m_accuracy * m_accuracy)
{
resched = true;
//last_resched_net = net;
}
resched_cnt++;
} while (resched && (resched_cnt < m_resched_loops / 3 ));
@ -146,7 +504,8 @@ ATTR_HOT inline int netlist_matrix_solver_t::solve_non_dynamic()
}
ATTR_HOT inline bool netlist_matrix_solver_t::solve()
template <int m_N, int _storage_N>
ATTR_HOT inline bool netlist_matrix_solver_gauss_seidel_t< m_N, _storage_N>::solve()
{
int resched_cnt = 0;
ATTR_UNUSED netlist_net_t *last_resched_net = NULL;
@ -171,17 +530,7 @@ ATTR_HOT inline bool netlist_matrix_solver_t::solve()
int this_resched;
do
{
/* update all non-linear devices */
for (netlist_core_device_t * const *p = m_dynamic.first(); p != NULL; p = m_dynamic.next(p))
switch ((*p)->family())
{
case netlist_device_t::DIODE:
static_cast<NETLIB_NAME(D) *>((*p))->update_terminals();
break;
default:
(*p)->update_terminals();
break;
}
update_dynamic();
this_resched = solve_non_dynamic();
resched_cnt += this_resched;
} while (this_resched > 1 && resched_cnt < m_resched_loops);
@ -194,32 +543,11 @@ ATTR_HOT inline bool netlist_matrix_solver_t::solve()
m_resched = true;
//if (resched)
//printf("Resched on net %s first term %s\n", last_resched_net->name().cstr(), last_resched_net->m_terms[0]->name().cstr());
//NL_VERBOSE_OUT(("Resched on net %s first term %s\n", last_resched_net->name().cstr(), last_resched_net->m_terms[0]->name().cstr());
return m_resched;
}
ATTR_HOT void netlist_matrix_solver_t::schedule()
{
if (!solve())
{
// printf("update_inputs\n");
update_inputs();
}
else
{
//printf("resched\n");
if (m_owner != NULL) this->m_owner->schedule1();
}
//solve();
// update_inputs();
}
ATTR_COLD void netlist_matrix_solver_t::reset()
{
m_last_step = netlist_time::zero;
}
// ----------------------------------------------------------------------------------------
// solver
// ----------------------------------------------------------------------------------------
@ -243,11 +571,14 @@ ATTR_COLD static void process_net(net_groups_t groups, int &cur_group, netlist_n
if (net->m_head == NULL)
return;
/* add the net */
NL_VERBOSE_OUT(("add %d - %s\n", cur_group, net->name().cstr()));
groups[cur_group].add(net);
for (netlist_core_terminal_t *p = net->m_head; p != NULL; p = p->m_update_list_next)
{
NL_VERBOSE_OUT(("terminal %s\n", p->name().cstr()));
if (p->isType(netlist_terminal_t::TERMINAL))
{
NL_VERBOSE_OUT(("isterminal\n"));
netlist_terminal_t *pt = static_cast<netlist_terminal_t *>(p);
netlist_net_t *other_net = &pt->m_otherterm->net();
if (!already_processed(groups, cur_group, other_net))
@ -327,12 +658,12 @@ NETLIB_UPDATE(solver)
{
NL_VERBOSE_OUT(("Step!\n"));
/* update all terminals for new time step */
m_last_step = now;
for (netlist_matrix_solver_t * const *e = m_mat_solvers.first(); e != NULL; e = m_mat_solvers.next(e))
{
(*e)->step(delta);
}
}
m_last_step = now;
#if HAS_OPENMP && USE_OPENMP
@ -388,6 +719,7 @@ ATTR_COLD void NETLIB_NAME(solver)::post_start()
// determine net groups
for (netlist_net_t * const *pn = netlist().m_nets.first(); pn != NULL; pn = netlist().m_nets.next(pn))
{
NL_VERBOSE_OUT(("proc %s\n", (*pn)->name().cstr()));
if (!already_processed(groups, cur_group, *pn))
{
cur_group++;
@ -399,7 +731,32 @@ ATTR_COLD void NETLIB_NAME(solver)::post_start()
SOLVER_VERBOSE_OUT(("Found %d net groups in %d nets\n", cur_group + 1, netlist().m_nets.count()));
for (int i = 0; i <= cur_group; i++)
{
netlist_matrix_solver_t *ms = new netlist_matrix_solver_t();
netlist_matrix_solver_t *ms;
switch (groups[i].count())
{
case 1:
ms = new netlist_matrix_solver_direct1_t();
break;
case 2:
//ms = new netlist_matrix_solver_gauss_seidel_t<2,2>();
// Below is only half as fast as gauss_seidel ....
ms = new netlist_matrix_solver_direct2_t();
break;
case 3:
ms = new netlist_matrix_solver_gauss_seidel_t<3,3>();
break;
case 4:
ms = new netlist_matrix_solver_gauss_seidel_t<4,4>();
break;
case 5:
ms = new netlist_matrix_solver_gauss_seidel_t<5,5>();
break;
default:
ms = new netlist_matrix_solver_gauss_seidel_t<0,16>();
break;
}
ms->m_accuracy = m_accuracy.Value();
ms->m_convergence_factor = m_convergence.Value();
ms->m_resched_loops = m_resched_loops.Value();

View File

@ -30,15 +30,18 @@ public:
typedef netlist_list_t<netlist_matrix_solver_t *> list_t;
typedef netlist_core_device_t::list_t dev_list_t;
netlist_matrix_solver_t() : m_resched(false), m_owner(NULL) {}
netlist_matrix_solver_t() : m_owner(NULL) {}
virtual ~netlist_matrix_solver_t() {}
ATTR_COLD void setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner);
ATTR_COLD virtual void setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner);
// return true if a reschedule is needed ...
ATTR_HOT bool solve();
ATTR_HOT int solve_non_dynamic();
ATTR_HOT void step(const netlist_time delta);
ATTR_HOT virtual bool solve() = 0;
ATTR_HOT virtual int solve_non_dynamic() = 0;
ATTR_HOT virtual void step(const netlist_time delta);
ATTR_HOT void update_inputs();
ATTR_HOT void update_dynamic();
ATTR_HOT void schedule();
@ -46,23 +49,123 @@ public:
ATTR_HOT inline bool is_timestep() { return m_steps.count() > 0; }
ATTR_HOT inline const NETLIB_NAME(solver) &owner() const;
ATTR_COLD void reset();
ATTR_COLD virtual void reset();
double m_accuracy;
double m_convergence_factor;
int m_resched_loops;
private:
protected:
netlist_net_t::list_t m_nets;
dev_list_t m_dynamic;
netlist_core_terminal_t::list_t m_inps;
dev_list_t m_steps;
bool m_resched;
netlist_time m_last_step;
NETLIB_NAME(solver) *m_owner;
};
template <int m_N, int _storage_N>
class netlist_matrix_solver_directb_t: public netlist_matrix_solver_t
{
public:
netlist_matrix_solver_directb_t() : netlist_matrix_solver_t() {}
virtual ~netlist_matrix_solver_directb_t() {}
ATTR_COLD virtual void setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner);
ATTR_HOT inline const int N() const { return (m_N == 0) ? m_nets.count() : m_N; }
ATTR_COLD virtual void reset() { netlist_matrix_solver_t::reset(); }
protected:
ATTR_HOT inline void build_LE(double (* RESTRICT m_A)[_storage_N], double (* RESTRICT m_RHS));
private:
ATTR_COLD int get_net_idx(netlist_net_t *net);
struct terms_t{
int net_this;
int net_other;
netlist_terminal_t *term;
};
int m_term_num;
int m_rail_start;
terms_t m_terms[100];
};
template <int m_N, int _storage_N>
class netlist_matrix_solver_gauss_seidel_t: public netlist_matrix_solver_t
{
public:
netlist_matrix_solver_gauss_seidel_t() : netlist_matrix_solver_t(), m_resched(false) {}
virtual ~netlist_matrix_solver_gauss_seidel_t() {}
ATTR_COLD virtual void setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner)
{
m_resched = true;
netlist_matrix_solver_t::setup(nets, owner);
}
// return true if a reschedule is needed ...
ATTR_HOT bool solve();
ATTR_HOT int solve_non_dynamic();
ATTR_HOT inline const int N() const { if (m_N == 0) return m_nets.count(); else return m_N; }
ATTR_COLD virtual void reset() { netlist_matrix_solver_t::reset(); }
private:
bool m_resched;
};
class netlist_matrix_solver_direct1_t: public netlist_matrix_solver_directb_t<1,1>
{
public:
netlist_matrix_solver_direct1_t() : netlist_matrix_solver_directb_t() {}
virtual ~netlist_matrix_solver_direct1_t() {}
ATTR_COLD virtual void setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner)
{
netlist_matrix_solver_directb_t::setup(nets, owner);
}
// return true if a reschedule is needed ...
ATTR_HOT bool solve();
ATTR_HOT int solve_non_dynamic();
//ATTR_COLD virtual void reset() { netlist_matrix_solver_t::reset(); }
private:
};
class netlist_matrix_solver_direct2_t: public netlist_matrix_solver_directb_t<2,2>
{
public:
netlist_matrix_solver_direct2_t() : netlist_matrix_solver_directb_t() {}
virtual ~netlist_matrix_solver_direct2_t() {}
ATTR_COLD virtual void setup(netlist_net_t::list_t &nets, NETLIB_NAME(solver) &owner)
{
netlist_matrix_solver_directb_t::setup(nets, owner);
}
// return true if a reschedule is needed ...
ATTR_HOT bool solve();
ATTR_HOT int solve_non_dynamic();
//ATTR_COLD virtual void reset() { netlist_matrix_solver_t::reset(); }
private:
};
NETLIB_DEVICE_WITH_PARAMS(solver,
typedef netlist_core_device_t::list_t dev_list_t;

View File

@ -33,8 +33,10 @@ NETLIB_UPDATE(twoterm)
netlist().solver()->schedule1();
#else
/* we only need to call the non-rail terminal */
m_P.net().solve();
m_N.net().solve();
if (!m_P.net().isRailNet())
m_P.net().solve();
else
m_N.net().solve();
#endif
}

View File

@ -86,6 +86,18 @@ public:
ATTR_HOT inline double deltaV() { return m_P.net().Q_Analog()- m_N.net().Q_Analog(); }
ATTR_HOT void set_mat(double a11, double a12, double a21, double a22, double r1, double r2)
{
m_P.m_gt = a11;
m_P.m_go = -a12;
m_N.m_gt = a22;
m_N.m_go = -a21;
m_P.m_Idr = -r1;
m_N.m_Idr = -r2;
}
protected:
ATTR_COLD virtual void start();
ATTR_COLD virtual void reset();
@ -183,19 +195,21 @@ public:
const double eVDVt = exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - 1.0);
m_G = m_Is * m_VtInv * eVDVt;
m_G = m_Is * m_VtInv * eVDVt + NETLIST_GMIN;
}
else
{
#if defined(_MSC_VER) && _MSC_VER < 1800
m_Vd = m_Vd + log((nVd - m_Vd) * m_VtInv + 1.0) * m_Vt;
#else
m_Vd = m_Vd + log1p((nVd - m_Vd) * m_VtInv) * m_Vt;
double a = (nVd - m_Vd) * m_VtInv;
if (a<1e-12 - 1.0) a = 1e-12 - 1.0;
m_Vd = m_Vd + log1p(a) * m_Vt;
#endif
const double eVDVt = exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - 1.0);
m_G = m_Is * m_VtInv * eVDVt;
m_G = m_Is * m_VtInv * eVDVt + NETLIST_GMIN;
}
//printf("nVd %f m_Vd %f Vcrit %f\n", nVd, m_Vd, m_Vcrit);
@ -212,9 +226,10 @@ public:
m_VtInv = 1.0 / m_Vt;
}
ATTR_HOT inline double I() { return m_Id; }
ATTR_HOT inline double G() { return m_G; }
ATTR_HOT inline double Ieq() { return (m_Id - m_Vd * m_G); }
ATTR_HOT inline double I() const { return m_Id; }
ATTR_HOT inline double G() const { return m_G; }
ATTR_HOT inline double Ieq() const { return (m_Id - m_Vd * m_G); }
ATTR_HOT inline double Vd() const { return m_Vd; }
/* owning object must save those ... */

View File

@ -553,7 +553,7 @@ public:
typedef netlist_list_t<netlist_net_t *> list_t;
friend class NETLIB_NAME(mainclock);
friend class netlist_matrix_solver_t;
friend class netlist_matrix_solver_t;
friend class netlist_logic_output_t;
friend class netlist_analog_output_t;
@ -638,6 +638,7 @@ public:
typedef netlist_list_t<netlist_terminal_t *> terminal_list_t;
terminal_list_t m_terms;
terminal_list_t m_rails; // FIXME: Make the solver use this !
netlist_matrix_solver_t *m_solver;
ATTR_HOT void solve();
@ -659,7 +660,7 @@ public:
netlist_object_t::save_register();
}
protected:
//protected: FIXME: needed by current solver code
UINT32 m_num_cons;
@ -667,6 +668,8 @@ protected:
hybrid_t m_cur;
hybrid_t m_new;
protected:
/* we don't use this to save state
* because we may get deleted again ...
*/

View File

@ -188,10 +188,10 @@ void netlist_setup_t::register_object(netlist_device_t &dev, const pstring &name
case netlist_param_t::LOGIC:
{
NL_VERBOSE_OUT(("Found parameter ... %s : %s\n", name.cstr(), val->cstr()));
int vald = 0;
if (sscanf(val.cstr(), "%d", &vald) != 1)
double vald = 0;
if (sscanf(val.cstr(), "%lf", &vald) != 1)
netlist().error("Invalid number conversion %s : %s\n", name.cstr(), val.cstr());
dynamic_cast<netlist_param_int_t &>(param).initial(vald);
dynamic_cast<netlist_param_int_t &>(param).initial((int) vald);
}
break;
case netlist_param_t::STRING: