Commit a7245160 authored by James Nutaro (1qn)'s avatar James Nutaro (1qn)
Browse files

Forgot the solver. Small code cleanup for m2c

parent a47f0541
......@@ -37,5 +37,5 @@ app: ${OBJS}
${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}
clean:
rm -f *.o a.out grid2ode func.cpp soln.txt m2c lex.yy.c *.gnu circuit.mo
rm -f *.o a.out grid2ode func.cpp soln.txt m2c lex.yy.c *.gnu circuit.mo sim vars.h
#include "sim.h"
#include <cfloat>
#include <algorithm>
#include <vector>
#include <list>
#include <cassert>
#include <iostream>
#include <cmath>
#include <cstring>
using namespace std;
// Current simulation time
static double tNow = 0.0;
// Maximum quanta of time or state
static const double step = 1E-6;
// Simulation ending time
static double tEnd = 0.0;
// Reporting interval
static double cint = 0.0;
// Smallest step size in reporting interval
static double hmin = 0.0;
// Model variable values
static double *q, *dq, *alg;
void print()
{
cout << tNow;
for (int i = 0; i < numStateVars; i++)
cout << " " << q[i];
for (int i = 0; i < numAlgVars; i++)
cout << " " << alg[i];
cout << " " << hmin << endl;
}
int main(int argc, char** argv)
{
int reports = 1;
double h, hReport, hEnd;
for (int i = 1; i < argc; i++)
{
if (strcmp("cint",argv[i]) == 0)
cint = atof(argv[++i]);
else if (strcmp("tend",argv[i]) == 0)
tEnd = atof(argv[++i]);
}
alg = new double[numAlgVars];
dq = new double[numStateVars];
q = new double[numStateVars];
// Calculate initial values
init_state(q);
event(dq,q,tNow);
func(dq,q,alg,0.0);
event(dq,q,tNow);
// Record the initial state
print();
hmin = DBL_MAX;
// Simulate until requested end of interval
do
{
// Get the smaller of the time to next report,
// and the stable time step. The step needs
// to be very small to resolve fast transients
// in the circuit
h = step;
for (int i = 0; i < numStateVars; i++)
if (dq[i] != 0.0)
h = ::min(h,step/fabs(dq[i]));
// We'll report the smallest step selected in the
// reporting interval
hmin = ::min(hmin,h);
hReport = cint*double(reports)-tNow;
hEnd = tEnd-tNow;
h = ::min(h,::min(hReport,hEnd));
// Update the state variables
for (int i = 0; i < numStateVars; i++)
q[i] += h*dq[i];
// Advance the clock
tNow += h;
// Apply events to state variable values
event(dq,q,tNow);
// Calculate new variable values for tNow
func(dq,q,alg,tNow);
// Apply derivatives values
event(dq,q,tNow);
// Report if its time or at the end of the simulation
if (h == hReport || h == hEnd)
{
reports++;
print();
hmin = DBL_MAX;
}
}
while (h != hEnd); // Stop at the end
// Cleanup and exit
delete [] dq;
delete [] q;
delete [] alg;
return 0;
}
......@@ -140,7 +140,6 @@ void number(const char* num)
%option yylineno
%option noyywrap
%option nodefault
ID [A-Za-z][A-Za-z_0-9]*
NUMBER -?[0-9]*\.?[0-9]+([e][-+]?[0-9]+)?
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment