Index: dc.c =================================================================== --- dc.c (revision 195) +++ dc.c (working copy) @@ -35,67 +35,23 @@ #include #include -/* assemble G,B,C,D,i,e into A and z */ -static void -CompileMatrices(ES_SimDC *sim) -{ - /* Compose the right-hand side vector (i, e). */ - M_Compose21(sim->z, sim->i, sim->e); - - /* Compose the coefficient matrix A from [G,B;C,D] and compute LU. */ - M_Compose22(sim->A, sim->G, sim->B, sim->C, sim->D); -} - -/* Formulate the KCL/branch equations and compute A=LU. */ +/* Solve Ax=z where A=[G,B;C,D], x=[v,j] and z=[i;e]. */ static int -LoadMatrices(ES_SimDC *sim, ES_Circuit *ckt) -{ - ES_Component *com; - - M_SetZero(sim->G); - M_SetZero(sim->B); - M_SetZero(sim->C); - M_SetZero(sim->D); - M_SetZero(sim->i); - M_SetZero(sim->e); - - /* Load datum node equation */ - - sim->G->v[0][0] = 1.0; - - /* Formulate the general equations. */ - CIRCUIT_FOREACH_COMPONENT(com, ckt) { - - if (com->loadDC_G != NULL) - com->loadDC_G(com, sim->G); - if (com->loadDC_BCD != NULL) - com->loadDC_BCD(com, sim->B, sim->C, sim->D); - if (com->loadDC_RHS != NULL) - com->loadDC_RHS(com, sim->i, sim->e); -#if 0 - if (com->loadAC != NULL) - com->loadAC(com, sim->Y); - if (com->loadSP != NULL) - com->loadSP(com, sim->S, sim->N); -#endif - } - - CompileMatrices(sim); - - return 0; -} - -static int SolveMNA(ES_SimDC *sim, ES_Circuit *ckt) { M_Real d; - M_FactorizeLU(sim->A, sim->LU, sim->piv, &d); + + M_Compose21(sim->z, sim->i, sim->e); + M_Compose22(sim->A, sim->G, sim->B, sim->C, sim->D); + /* Find LU factorization and solve by backsubstitution. */ + M_FactorizeLU(sim->A, sim->LU, sim->piv, &d); M_Copy(sim->x, sim->z); M_BacksubstLU(sim->LU, sim->piv, sim->x); return (0); } +/* Simulation timestep. */ static Uint32 StepMNA(void *obj, Uint32 ival, void *arg) { @@ -103,10 +59,12 @@ ES_SimDC *sim = arg; ES_Component *com; static Uint32 t1 = 0; - M_Vector *xprev; + M_Vector *xPrev; M_Real diff; Uint i = 1, j; + xPrev = M_New(sim->x->m, 1); + if (ckt->m == 0) { AG_SetError(_("Circuit has no voltage sources")); goto halt; @@ -118,70 +76,62 @@ AG_PostEvent(NULL, ckt, "circuit-step-begin", NULL); - /* Invoke the time step function of all models. */ + /* Formulate and solve for the initial conditions. */ CIRCUIT_FOREACH_COMPONENT(com, ckt) { - if (com->intStep != NULL) - com->intStep(com, sim->Telapsed); - if (com->dcStep != NULL) - com->dcUpdate(com, sim->G, sim->B, sim->C, sim->D, sim->i, sim->e, sim->v, sim->j); + if (com->dcStepBegin != NULL) + com->dcStepBegin(com, sim->G, sim->B, sim->C, sim->D, + sim->i, sim->e, sim->v, sim->j); } - CompileMatrices(sim); - if (SolveMNA(sim, ckt) == -1) - { goto halt; - } - do - { - if (++i > sim->max_iters) { + /* Iterate until a stable solution is found. */ + do { + if (++i > sim->itersMax) { AG_SetError(_("Could not find stable solution in " "%u iterations"), i); goto halt; } - - /* - * Save the A matrix and execute per-component procedural update - * functions. The intUpdate function can modify either A or the - * RHS in response to calculated node voltages. - */ - xprev = M_Dup(sim->x); + M_Copy(xPrev, sim->x); CIRCUIT_FOREACH_COMPONENT(com, ckt) { - if (com->intUpdate != NULL) - com->intUpdate(com); - if (com->dcUpdate != NULL) - com->dcUpdate(com, sim->G, sim->B, sim->C, sim->D, sim->i, sim->e, sim->v, sim->j); + if (com->dcStepIter != NULL) + com->dcStepIter(com, sim->G, sim->B, sim->C, + sim->D, sim->i, sim->e, sim->v, sim->j); } - CompileMatrices(sim); - if (SolveMNA(sim, ckt) == -1) - { goto halt; - } /* compute error */ - diff = 0; - for (j=0; j < sim->x->m; j++) - { - diff += (sim->x->v[j][0] - xprev->v[j][0])*(sim->x->v[j][0] - xprev->v[j][0]); + for (j = 0; j < sim->x->m; j++) { + diff += (sim->x->v[j][0] - xPrev->v[j][0]) * + (sim->x->v[j][0] - xPrev->v[j][0]); } - diff = sqrt(diff); - + diff = Sqrt(diff); } while (diff > 0.001); + + /* Invoke the Component DC specific post-timestep callbacks. */ + CIRCUIT_FOREACH_COMPONENT(com, ckt) { + if (com->dcStepEnd != NULL) + com->dcStepEnd(com, sim->A, sim->x, sim->z); + } - if (i > sim->iters_hiwat) { sim->iters_hiwat = i; } - else if (i < sim->iters_lowat) { sim->iters_lowat = i; } - + /* Update the statistics. */ + if (i > sim->itersHiwat) { sim->itersLowat = i; } sim->TavgReal = SDL_GetTicks() - t1; t1 = SDL_GetTicks(); sim->Telapsed++; + + /* Invoke the general post-timestep callbacks (used by Scope, etc.) */ AG_PostEvent(NULL, ckt, "circuit-step-end", NULL); + + M_Free(xPrev); return (ival); halt: AG_TextMsg(AG_MSG_ERROR, _("%s; simulation stopped"), AG_GetError()); - SIM(sim)->running = 0; + StopSimulation(sim); + M_Free(xPrev); return (0); } @@ -194,11 +144,13 @@ sim->Telapsed = 0; sim->TavgReal = 0; sim->speed = 500; - sim->max_iters = 10000; - sim->iters_hiwat = 1; - sim->iters_lowat = 1; - AG_SetTimeout(&sim->update_to, StepMNA, sim, AG_CANCEL_ONDETACH); + sim->itersMax = 10000; + sim->itersHiwat = 1; + sim->itersLowat = 1; + AG_SetTimeout(&sim->toUpdate, StepMNA, sim, AG_CANCEL_ONDETACH); + sim->T0 = 290.0; + sim->A = M_New(0,0); sim->G = M_New(0,0); sim->B = M_New(0,0); @@ -260,35 +212,55 @@ { ES_SimDC *sim = p; ES_Circuit *ckt = SIM(sim)->ckt; - + ES_Component *com; + + /* Reallocate the matrices if necessary. */ Resize(sim, ckt); - + + /* Register the real-time simulation timer. */ SIM(sim)->running = 1; AG_LockTimeouts(ckt); - if (AG_TimeoutIsScheduled(ckt, &sim->update_to)) { - AG_DelTimeout(ckt, &sim->update_to); + if (AG_TimeoutIsScheduled(ckt, &sim->toUpdate)) { + AG_DelTimeout(ckt, &sim->toUpdate); } - AG_AddTimeout(ckt, &sim->update_to, 1000/sim->speed); + AG_AddTimeout(ckt, &sim->toUpdate, 1000/sim->speed); AG_UnlockTimeouts(ckt); sim->Telapsed = 0; - LoadMatrices(sim, ckt); + /* Initialize the matrices. */ + M_SetZero(sim->G); + M_SetZero(sim->B); + M_SetZero(sim->C); + M_SetZero(sim->D); + M_SetZero(sim->i); + M_SetZero(sim->e); + /* Load datum node equation */ + sim->G->v[0][0] = 1.0; + + AG_PostEvent(NULL, ckt, "circuit-sim-begin", NULL); ES_SimLog(sim, _("Simulation started (m=%u, n=%u)"), ckt->m, ckt->n); + return 0; } static void +StopSimulation(ES_SimDC *sim) +{ + SIM(sim)->running = 0; + ES_SimLog(sim, _("Simulation stopped at %uns."), sim->Telapsed); + AG_PostEvent(NULL, ckt, "circuit-sim-end", NULL); +} + +static void Stop(void *p) { ES_SimDC *sim = p; ES_Circuit *ckt = SIM(sim)->ckt; AG_LockTimeouts(ckt); - AG_DelTimeout(ckt, &sim->update_to); + AG_DelTimeout(ckt, &sim->toUpdate); AG_UnlockTimeouts(ckt); - SIM(sim)->running = 0; - - ES_SimLog(sim, _("Simulation stopped at %uns."), sim->Telapsed); + StopSimulation(sim); } static void @@ -355,7 +327,7 @@ sbu = AG_SpinbuttonNew(ntab, 0, _("Maximum iterations per step: ")); AG_WidgetBind(sbu, "value", AG_WIDGET_UINT, - &sim->max_iters); + &sim->itersMax); sbu = AG_SpinbuttonNew(ntab, 0, _("Speed (updates/sec): ")); AG_WidgetBind(sbu, "value", AG_WIDGET_UINT32, &sim->speed); @@ -369,7 +341,7 @@ &sim->TavgReal); AG_LabelNewPolled(ntab, 0, _("Iterations/step watermark: %u-%u"), - &sim->iters_lowat, &sim->iters_hiwat); + &sim->itersLowat, &sim->itersHiwat); } ntab = AG_NotebookAddTab(nb, "[A]", AG_BOX_VERT);