Line data Source code
1 : // Copyright 2014 Bobby Powers. All rights reserved.
2 : // Use of this source code is governed by a BSD-style
3 : // license that can be found in the LICENSE file.
4 :
5 : #include <math.h>
6 : #include <stdlib.h>
7 : #include <string.h>
8 :
9 : #include "utf.h"
10 : #include "sd.h"
11 : #include "sd_internal.h"
12 :
13 :
14 : typedef struct {
15 : Walker w;
16 : AVar *module;
17 : AVar *av;
18 : Node *curr;
19 : } AVarWalker;
20 :
21 : typedef struct {
22 : const char *const name;
23 : Fn fn;
24 : } FnDef;
25 :
26 : static double rt_min(SDSim *s, Node *n, double dt, double t, size_t len, double *args);
27 : static double rt_max(SDSim *s, Node *n, double dt, double t, size_t len, double *args);
28 : static double rt_pulse(SDSim *s, Node *n, double dt, double t, size_t len, double *args);
29 :
30 : static double *sim_curr(SDSim *s);
31 : static double *sim_next(SDSim *s);
32 :
33 : static void calc(SDSim *s, double *data, Slice *l, bool initial);
34 : static void calc_stocks(SDSim *s, double *data, Slice *l);
35 :
36 : static double svisit(SDSim *s, Node *n, double dt, double time);
37 :
38 : static AVar *module(SDProject *p, AVar *parent, SDModel *model, Var *module);
39 : static int module_compile(AVar *module);
40 : static int module_assign_offsets(AVar *module, int *offset);
41 : static int module_get_varnames(AVar *module, const char **result, size_t max);
42 : static void module_clear_visited(AVar *module);
43 : static int module_sort_runlists(AVar *module);
44 : static int module_add_to_runlists(AVar *module, AVar *av);
45 :
46 : static const char *avar_qual_name(AVar *av);
47 :
48 : static AVarWalker *avar_walker_new(AVar *module, AVar *av);
49 : static void avar_walker_ref(void *data);
50 : static void avar_walker_unref(void *data);
51 : static void avar_walker_start(void *data, Node *n);
52 : static Walker *avar_walker_start_child(void *data, Node *n);
53 : static void avar_walker_end_child(void *data, Node *);
54 :
55 : static const WalkerOps AVAR_WALKER_OPS = {
56 : .ref = avar_walker_ref,
57 : .unref = avar_walker_unref,
58 : .start = avar_walker_start,
59 : .start_child = avar_walker_start_child,
60 : .end_child = avar_walker_end_child,
61 : .end = NULL,
62 : };
63 :
64 : static const FnDef RT_FNS[] = {
65 : {"pulse", rt_pulse},
66 : {"min", rt_min},
67 : {"max", rt_max},
68 : };
69 : static const size_t RT_FNS_LEN = sizeof(RT_FNS)/sizeof(RT_FNS[0]);
70 :
71 : AVar *
72 100 : avar(AVar *parent, Var *v)
73 : {
74 : SDModel *model;
75 : AVar *av;
76 : int err;
77 :
78 100 : av = NULL;
79 100 : err = 0;
80 :
81 100 : if (v->type == VAR_MODULE) {
82 4 : model = sd_project_get_model(parent->model->file->project, v->name);
83 4 : if (!model)
84 0 : goto error;
85 4 : av = module(parent->model->file->project, parent, model, v);
86 4 : if (av)
87 4 : av->v = v;
88 4 : return av;
89 : }
90 :
91 96 : av = calloc(1, sizeof(*av));
92 96 : if (!av)
93 0 : goto error;
94 96 : av->v = v;
95 96 : av->parent = parent;
96 :
97 96 : if (v->eqn)
98 88 : err = avar_eqn_parse(av);
99 :
100 96 : if (err)
101 0 : goto error;
102 96 : av->is_const = av->node && av->node->type == N_FLOATLIT;
103 :
104 96 : return av;
105 : error:
106 0 : printf("eqn parse failed for %s\n", v->name);
107 0 : free(av);
108 0 : return NULL;
109 : }
110 :
111 : const char *
112 85 : avar_qual_name(AVar *av)
113 : {
114 : const char *parent;
115 : size_t parent_len, name_len;
116 :
117 85 : if (av->qual_name)
118 0 : return av->qual_name;
119 :
120 : // if we _are_ the root model, return <main>
121 85 : if (!av->parent)
122 0 : return "<main>";
123 :
124 : // if we're in the root model, return the var's name
125 85 : if (av->parent->parent == NULL)
126 73 : return av->v->name;
127 :
128 12 : parent = avar_qual_name(av->parent);
129 12 : if (!parent)
130 : // TODO(bp) handle ENOMEM
131 0 : return NULL;
132 12 : parent_len = strlen(parent);
133 :
134 12 : name_len = strlen(av->v->name);
135 :
136 12 : av->qual_name = malloc(parent_len + name_len + 2);
137 12 : if (!av->qual_name)
138 : // TODO(bp) handle ENOMEM
139 0 : return NULL;
140 :
141 12 : memcpy(av->qual_name, parent, parent_len);
142 12 : av->qual_name[parent_len] = '.';
143 12 : memcpy(&av->qual_name[parent_len+1], av->v->name, name_len);
144 12 : av->qual_name[parent_len+name_len+1] = '\0';
145 :
146 12 : return av->qual_name;
147 : }
148 :
149 : AVarWalker *
150 88 : avar_walker_new(AVar *module, AVar *av)
151 : {
152 88 : AVarWalker *w = calloc(1, sizeof(*w));
153 :
154 88 : if (!w)
155 0 : return NULL;
156 :
157 88 : w->w.ops = &AVAR_WALKER_OPS;
158 88 : avar_walker_ref(w);
159 88 : w->av = av;
160 88 : w->module = module;
161 :
162 88 : return w;
163 : }
164 :
165 : void
166 292 : avar_walker_ref(void *data)
167 : {
168 292 : AVarWalker *w = data;
169 292 : __sync_fetch_and_add(&w->w.refcount, 1);
170 292 : }
171 :
172 : void
173 292 : avar_walker_unref(void *data)
174 : {
175 292 : AVarWalker *w = data;
176 :
177 292 : if (!w || __sync_sub_and_fetch(&w->w.refcount, 1) != 0)
178 204 : return;
179 :
180 88 : free(w);
181 : }
182 :
183 : void
184 204 : avar_walker_start(void *data, Node *n)
185 : {
186 204 : AVarWalker *avw = data;
187 : AVar *dep;
188 :
189 204 : avw->curr = n;
190 :
191 204 : switch (n->type) {
192 : case N_IDENT:
193 75 : dep = resolve(avw->module, n->sval);
194 75 : if (!dep) {
195 0 : printf("resolve failed for %s\n", n->sval);
196 : // TODO: error
197 0 : return;
198 : }
199 75 : n->av = dep;
200 75 : slice_append(&avw->av->direct_deps, dep);
201 75 : break;
202 : case N_FLOATLIT:
203 72 : n->fval = atof(n->sval);
204 72 : break;
205 : case N_CALL:
206 7 : for (size_t i = 0; i < RT_FNS_LEN; i++) {
207 7 : if (strcmp(n->left->sval, RT_FNS[i].name) == 0) {
208 4 : n->fn = RT_FNS[i].fn;
209 4 : break;
210 : }
211 : }
212 4 : if (!n->fn)
213 0 : printf("unknown fn '%s' for call\n", n->left->sval);
214 4 : break;
215 : default:
216 53 : break;
217 : }
218 : }
219 :
220 : Walker *
221 120 : avar_walker_start_child(void *data, Node *n)
222 : {
223 120 : AVarWalker *avw = data;
224 120 : Node *curr = avw->curr;
225 :
226 : // skip trying to resolve function name identifiers for calls
227 : // a second time - already handled in avar_walker_start above
228 120 : if (curr->type == N_CALL && n == curr->left)
229 4 : return NULL;
230 :
231 116 : avw->w.ops->ref(&avw->w);
232 116 : return &avw->w;
233 : }
234 :
235 : void
236 116 : avar_walker_end_child(void *data, Node *n)
237 : {
238 116 : }
239 :
240 : int
241 132 : avar_init(AVar *av, AVar *module)
242 : {
243 : AVarWalker *w;
244 : bool ok;
245 :
246 132 : w = NULL;
247 :
248 : // is amodule if we have a model pointer
249 132 : if (av->model) {
250 20 : return module_compile(av);
251 112 : } else if (av->v->type == VAR_REF) {
252 8 : AVar *src = resolve(module->parent, av->v->src);
253 8 : if (src) {
254 8 : av->src = src;
255 8 : return 0;
256 : }
257 0 : goto error;
258 104 : } else if (!av->v->eqn && av->v->name && strcmp(av->v->name, "time") == 0) {
259 16 : return 0;
260 : }
261 :
262 88 : w = avar_walker_new(module, av);
263 :
264 88 : ok = node_walk(&w->w, av->node);
265 88 : if (!ok)
266 0 : goto error;
267 :
268 98 : for (size_t i = 0; i < av->v->inflows.len; i++) {
269 10 : char *in_name = av->v->inflows.elems[i];
270 10 : AVar *in = resolve(module, in_name);
271 10 : if (!in)
272 0 : goto error;
273 10 : slice_append(&av->inflows, in);
274 : }
275 101 : for (size_t i = 0; i < av->v->outflows.len; i++) {
276 13 : char *out_name = av->v->outflows.elems[i];
277 13 : AVar *out = resolve(module, out_name);
278 13 : if (!out)
279 0 : goto error;
280 13 : slice_append(&av->outflows, out);
281 : }
282 :
283 88 : avar_walker_unref(w);
284 :
285 88 : return SD_ERR_NO_ERROR;
286 : error:
287 0 : avar_walker_unref(w);
288 0 : return SD_ERR_UNSPECIFIED;
289 : }
290 :
291 : void
292 133 : avar_free(AVar *av)
293 : {
294 133 : if (!av)
295 1 : return;
296 :
297 132 : sd_model_unref(av->model);
298 248 : for (size_t i = 0; i < av->avars.len; i++) {
299 116 : AVar *child = av->avars.elems[i];
300 116 : avar_free(child);
301 : }
302 132 : free(av->avars.elems);
303 132 : free(av->qual_name);
304 132 : node_free(av->node);
305 132 : free(av->direct_deps.elems);
306 132 : free(av->inflows.elems);
307 132 : free(av->outflows.elems);
308 132 : free(av->initials.elems);
309 132 : free(av->flows.elems);
310 132 : free(av->stocks.elems);
311 132 : var_free(av->time);
312 132 : free(av);
313 : }
314 :
315 : AVar *
316 20 : module(SDProject *p, AVar *parent, SDModel *model, Var *vmodule)
317 : {
318 : AVar *module;
319 : Slice conns;
320 :
321 20 : module = calloc(1, sizeof(*module));
322 20 : if (!module)
323 0 : goto error;
324 20 : module->v = vmodule;
325 20 : module->model = model;
326 20 : module->parent = parent;
327 :
328 20 : if (!parent) {
329 16 : Var *time = calloc(1, sizeof(*time));
330 16 : if (!time)
331 0 : goto error;
332 16 : time->type = VAR_AUX;
333 16 : time->name = strdup("time");
334 16 : module->time = time;
335 :
336 16 : AVar *atime = calloc(1, sizeof(*atime));
337 16 : if (!atime)
338 0 : goto error;
339 16 : atime->v = time;
340 16 : atime->parent = module;
341 16 : slice_append(&module->avars, atime);
342 : }
343 :
344 20 : memset(&conns, 0, sizeof(conns));
345 20 : if (vmodule)
346 4 : conns = vmodule->conns;
347 :
348 120 : for (size_t i = 0; i < model->vars.len; i++) {
349 : AVar *av;
350 : Var *v;
351 :
352 100 : av = NULL;
353 100 : v = model->vars.elems[i];
354 : // FIXME: this is unelegant and seems inefficient.
355 152 : for (size_t j = 0; j < conns.len; j++) {
356 60 : Var *r = conns.elems[j];
357 60 : if (strcmp(v->name, r->name) == 0) {
358 8 : av = avar(module, r);
359 8 : break;
360 : }
361 : }
362 100 : if (!av)
363 92 : av = avar(module, v);
364 100 : if (!av)
365 0 : goto error;
366 100 : slice_append(&module->avars, av);
367 : }
368 :
369 20 : return module;
370 : error:
371 0 : avar_free(module);
372 0 : return NULL;
373 : }
374 :
375 : AVar *
376 189 : resolve(AVar *module, const char *name)
377 : {
378 : size_t len;
379 : const char *subvar;
380 :
381 189 : len = 0;
382 :
383 189 : if (!name)
384 0 : return NULL;
385 :
386 : // a name like .area gets resolved to area
387 189 : if (name[0] == '.')
388 4 : name++;
389 :
390 189 : subvar = strchr(name, '.');
391 189 : if (subvar) {
392 17 : len = subvar - name;
393 17 : subvar++;
394 : }
395 :
396 787 : for (size_t i = 0; i < module->avars.len; i++) {
397 785 : AVar *av = module->avars.elems[i];
398 785 : if (subvar && av->v->type == VAR_MODULE && strncmp(av->v->name, name, len) == 0) {
399 17 : return resolve(av, subvar);
400 : }
401 768 : else if (strcmp(av->v->name, name) == 0)
402 170 : return av;
403 : }
404 :
405 2 : return NULL;
406 : }
407 :
408 : int
409 20 : module_compile(AVar *module)
410 : {
411 : AVar *av;
412 : int err, failed;
413 :
414 20 : failed = 0;
415 136 : for (size_t i = 0; i < module->avars.len; i++) {
416 116 : av = module->avars.elems[i];
417 116 : err = avar_init(av, module);
418 116 : if (err)
419 0 : failed++;
420 : }
421 20 : if (failed)
422 0 : return SD_ERR_UNSPECIFIED;
423 :
424 136 : for (size_t i = 0; i < module->avars.len; i++) {
425 116 : av = module->avars.elems[i];
426 : // TODO: fix for multiple levels of indirection.
427 116 : if (av->v->type == VAR_REF)
428 8 : av->offset = av->src->offset;
429 : }
430 :
431 : // sorting of runlists is done in a separate step after
432 : // offsets are assigned.
433 :
434 20 : return SD_ERR_NO_ERROR;
435 : }
436 :
437 : SDSim *
438 17 : sd_sim_new(SDProject *p, const char *model_name)
439 : {
440 : SDSim *sim;
441 : SDModel *model;
442 : int err, offset;
443 :
444 17 : offset = 0;
445 17 : model = NULL;
446 17 : sim = calloc(1, sizeof(*sim));
447 17 : if (!sim)
448 0 : goto error;
449 17 : sd_sim_ref(sim);
450 :
451 : // FIXME: check refcounting
452 17 : model = sd_project_get_model(p, model_name);
453 17 : if (!model)
454 1 : goto error;
455 :
456 16 : sim->module = module(p, NULL, model, NULL);
457 16 : if (!sim->module)
458 0 : goto error;
459 :
460 16 : sd_project_ref(p);
461 16 : sim->project = p;
462 :
463 16 : err = avar_init(sim->module, NULL);
464 16 : if (err)
465 0 : goto error;
466 :
467 16 : err = module_assign_offsets(sim->module, &offset);
468 16 : if (err)
469 0 : goto error;
470 :
471 16 : err = module_sort_runlists(sim->module);
472 16 : if (err)
473 0 : goto error;
474 :
475 16 : sim->nvars = offset;
476 16 : err = sd_sim_reset(sim);
477 16 : if (err)
478 0 : goto error;
479 :
480 16 : return sim;
481 : error:
482 1 : sd_sim_unref(sim);
483 1 : return NULL;
484 : }
485 :
486 : int
487 20 : module_assign_offsets(AVar *module, int *offset)
488 : {
489 : int err;
490 136 : for (size_t i = 0; i < module->avars.len; i++) {
491 : AVar *av;
492 116 : av = module->avars.elems[i];
493 116 : if (av->model) {
494 4 : err = module_assign_offsets(av, offset);
495 4 : if (err)
496 0 : return err;
497 112 : } else if (!av->src) {
498 : // assign offsets for everything thats not a
499 : // module or ref
500 104 : av->offset = (*offset)++;
501 : }
502 : }
503 :
504 20 : return 0;
505 : }
506 :
507 106 : int module_add_to_runlists(AVar *module, AVar *av)
508 : {
509 106 : if (av->visited)
510 0 : return 0;
511 :
512 : // TODO: better circularity error reporting
513 106 : if (av->visiting)
514 0 : return SD_ERR_CIRCULAR;
515 :
516 106 : av->visiting = true;
517 :
518 : // make sure any of our dependencies are on runlists before
519 : // us.
520 181 : for (size_t i = 0; i < av->direct_deps.len; i++) {
521 75 : AVar *dep = av->direct_deps.elems[i];
522 : int err;
523 :
524 75 : if (dep->visited)
525 28 : continue;
526 :
527 47 : err = module_add_to_runlists(module, dep);
528 47 : if (err)
529 0 : return err;
530 : }
531 :
532 106 : if (av->v->type == VAR_MODULE) {
533 4 : slice_append(&module->initials, av);
534 4 : slice_append(&module->flows, av);
535 4 : slice_append(&module->stocks, av);
536 102 : } else if (av->v->type == VAR_STOCK) {
537 13 : slice_append(&module->initials, av);
538 13 : slice_append(&module->stocks, av);
539 : // refs are not simulated
540 89 : } else if (av->v->type != VAR_REF) {
541 81 : slice_append(&module->initials, av);
542 81 : if (av->is_const)
543 27 : slice_append(&module->stocks, av);
544 : else
545 54 : slice_append(&module->flows, av);
546 : }
547 :
548 106 : av->visited = true;
549 106 : av->visiting = false;
550 :
551 106 : return 0;
552 : }
553 :
554 : int
555 20 : module_sort_runlists(AVar *module)
556 : {
557 : int off;
558 :
559 20 : module_clear_visited(module);
560 :
561 20 : module->visiting = true;
562 :
563 20 : if (!module->parent)
564 16 : off = 1;
565 : else
566 4 : off = 0;
567 :
568 120 : for (size_t i = off; i < module->avars.len; i++) {
569 100 : AVar *sub = module->avars.elems[i];
570 : int err;
571 :
572 100 : if (sub->visited)
573 41 : continue;
574 :
575 59 : if (sub->v->type == VAR_MODULE) {
576 4 : err = module_sort_runlists(sub);
577 4 : if (err)
578 0 : return err;
579 : }
580 :
581 59 : err = module_add_to_runlists(module, sub);
582 59 : if (err)
583 0 : return err;
584 : }
585 :
586 20 : module->visiting = false;
587 :
588 20 : return 0;
589 : }
590 :
591 : int
592 17 : sd_sim_reset(SDSim *s)
593 : {
594 17 : int err = 0;
595 : size_t save_every;
596 :
597 17 : s->spec = s->module->model->file->sim_specs;
598 17 : s->step = 0;
599 17 : s->save_step = 0;
600 17 : s->nsteps = (s->spec.stop - s->spec.start)/s->spec.dt + 1;
601 :
602 17 : save_every = s->spec.savestep/s->spec.dt+.5;
603 17 : s->save_every = save_every > 1 ? save_every : 1;
604 17 : s->nsaves = s->nsteps/s->save_every;
605 17 : if (s->nsteps % s->save_every)
606 2 : s->nsaves++;
607 :
608 17 : free(s->slab);
609 : // XXX: 1 extra step to simplify run_to
610 17 : s->slab = calloc(s->nvars*(s->nsaves + 1), sizeof(double));
611 17 : if (!s->slab) {
612 0 : err = SD_ERR_NOMEM;
613 0 : goto error;
614 : }
615 17 : s->curr = s->slab;
616 17 : s->next = NULL;
617 :
618 17 : s->curr[TIME] = s->spec.start;
619 :
620 17 : calc(s, s->curr, &s->module->initials, true);
621 : error:
622 17 : return err;
623 : }
624 :
625 : void
626 1004878 : calc(SDSim *s, double *data, Slice *l, bool initial)
627 : {
628 : double dt;
629 :
630 1004878 : dt = s->spec.dt;
631 :
632 : //printf("CALC\n");
633 1016985 : for (size_t i = 0; i < l->len; i++) {
634 12107 : AVar *av = l->elems[i];
635 12107 : if (!av->node) {
636 374 : if (initial)
637 10 : calc(s, data, &av->initials, true);
638 : else
639 364 : calc(s, data, &av->flows, false);
640 374 : continue;
641 : }
642 11733 : double v = svisit(s, av->node, dt, data[0]);
643 11733 : if (av->v->gf)
644 1728 : v = lookup(av->v->gf, v);
645 11733 : data[av->offset] = v;
646 : }
647 1004878 : }
648 :
649 : void
650 1004583 : calc_stocks(SDSim *s, double *data, Slice *l)
651 : {
652 : double prev, v, dt;
653 :
654 1004583 : dt = s->spec.dt;
655 :
656 : //printf("CALC STOCKS\n");
657 4029180 : for (size_t i = 0; i < l->len; i++) {
658 3024597 : AVar *av = l->elems[i];
659 : // XXX: this could also be implemented by building the
660 : // addition and subtraction of flows in the stock's
661 : // AST. Maybe that would be cleaner?
662 3024597 : switch (av->v->type) {
663 : case VAR_STOCK:
664 1011325 : prev = s->curr[av->offset];
665 1011325 : v = 0;
666 2018967 : for (size_t i = 0; i < av->inflows.len; i++) {
667 1007642 : AVar *in = av->inflows.elems[i];
668 1007642 : v += s->curr[in->offset];
669 : }
670 1019796 : for (size_t i = 0; i < av->outflows.len; i++) {
671 8471 : AVar *out = av->outflows.elems[i];
672 8471 : v -= s->curr[out->offset];
673 : }
674 1011325 : data[av->offset] = prev + v*dt;
675 1011325 : break;
676 : case VAR_MODULE:
677 96 : calc_stocks(s, data, &av->stocks);
678 96 : break;
679 : default:
680 2013176 : v = svisit(s, av->node, dt, s->curr[0]);
681 2013176 : data[av->offset] = v;
682 2013176 : break;
683 : }
684 : }
685 1004583 : }
686 :
687 : int
688 17 : sd_sim_run_to(SDSim *s, double end)
689 : {
690 : double dt;
691 :
692 17 : if (!s)
693 1 : return -1;
694 :
695 16 : dt = s->spec.dt;
696 16 : s->curr = sim_curr(s);
697 16 : s->next = sim_next(s);
698 :
699 1004503 : while (s->step < s->nsteps && s->curr[TIME] <= end) {
700 1004487 : calc(s, s->curr, &s->module->flows, false);
701 1004487 : calc_stocks(s, s->next, &s->module->stocks);
702 :
703 1004487 : if (s->step + 1 == s->nsteps)
704 16 : break;
705 :
706 : // calculate this way instead of += dt to minimize
707 : // cumulative floating point errors.
708 1004471 : s->next[TIME] = s->spec.start + (s->step+1)*dt;
709 :
710 1004471 : if (s->step++ % s->save_every != 0) {
711 999990 : memcpy(s->curr, s->next, s->nvars*sizeof(double));
712 : } else {
713 4481 : s->save_step++;
714 4481 : s->curr = sim_curr(s);
715 4481 : s->next = sim_next(s);
716 : }
717 : }
718 :
719 16 : return 0;
720 : }
721 :
722 : int
723 17 : sd_sim_run_to_end(SDSim *s)
724 : {
725 17 : if (!s)
726 1 : return -1;
727 16 : return sd_sim_run_to(s, s->spec.stop + 1);
728 : }
729 :
730 : double *
731 4497 : sim_curr(SDSim *s)
732 : {
733 4497 : return &s->slab[s->save_step*s->nvars];
734 : }
735 :
736 : double *
737 4497 : sim_next(SDSim *s)
738 : {
739 4497 : return &s->slab[(s->save_step+1)*s->nvars];
740 : }
741 :
742 : void
743 18 : sd_sim_ref(SDSim *sim)
744 : {
745 18 : if (!sim)
746 1 : return;
747 17 : __sync_fetch_and_add(&sim->refcount, 1);
748 : }
749 :
750 : int
751 7 : sd_sim_get_value(SDSim *s, const char *name, double *result)
752 : {
753 : AVar *av;
754 :
755 7 : if (!s || !name || !result)
756 1 : return SD_ERR_UNSPECIFIED;
757 :
758 6 : if (strcmp(name, "time") == 0) {
759 1 : *result = s->curr[TIME];
760 1 : return 0;
761 : }
762 :
763 5 : av = resolve(s->module, name);
764 5 : if (!av)
765 1 : return SD_ERR_UNSPECIFIED;
766 :
767 4 : *result = s->curr[av->offset];
768 4 : return 0;
769 : }
770 :
771 : void
772 17 : sd_sim_unref(SDSim *sim)
773 : {
774 17 : if (!sim)
775 0 : return;
776 17 : if (__sync_sub_and_fetch(&sim->refcount, 1) == 0) {
777 17 : avar_free(sim->module);
778 17 : sd_project_unref(sim->project);
779 17 : free(sim->slab);
780 17 : free(sim);
781 : }
782 : }
783 :
784 : double
785 2058940 : svisit(SDSim *s, Node *n, double dt, double time)
786 : {
787 2058940 : double v = NAN;
788 : double cond, l, r;
789 : double args[6];
790 : int off;
791 :
792 2058940 : switch (n->type) {
793 : case N_PAREN:
794 486 : v = svisit(s, n->left, dt, time);
795 486 : break;
796 : case N_FLOATLIT:
797 2013567 : v = n->fval;
798 2013567 : break;
799 : case N_IDENT:
800 28132 : if (n->av->src)
801 200 : off = n->av->src->offset;
802 : else
803 27932 : off = n->av->offset;
804 28132 : v = s->curr[off];
805 28132 : break;
806 : case N_CALL:
807 74 : memset(args, 0, 6*sizeof(*args));
808 74 : (void)n->left->sval;
809 272 : for (size_t i = 0; i < n->args.len; i++) {
810 198 : Node *arg = n->args.elems[i];
811 198 : args[i] = svisit(s, arg, dt, time);
812 : }
813 74 : v = n->fn(s, n, dt, time, n->args.len, args);
814 74 : break;
815 : case N_IF:
816 68 : cond = svisit(s, n->cond, dt, time);
817 68 : if (cond != 0)
818 43 : v = svisit(s, n->left, dt, time);
819 : else
820 25 : v = svisit(s, n->right, dt, time);
821 68 : break;
822 : case N_UNARY:
823 15 : l = svisit(s, n->left, dt, time);
824 15 : switch (n->op) {
825 : case '+':
826 0 : v = l;
827 0 : break;
828 : case '-':
829 12 : v = -l;
830 12 : break;
831 : case '!':
832 3 : v = l == 0 ? 1 : 0;
833 3 : break;
834 : }
835 15 : break;
836 : case N_BINARY:
837 16598 : l = svisit(s, n->left, dt, time);
838 16598 : r = svisit(s, n->right, dt, time);
839 16598 : switch (n->op) {
840 : case '+':
841 5 : v = l + r;
842 5 : break;
843 : case '-':
844 489 : v = l - r;
845 489 : break;
846 : case '*':
847 8532 : v = l * r;
848 8532 : break;
849 : case '/':
850 7423 : v = l / r;
851 7423 : break;
852 : case '<':
853 12 : v = l < r ? 1 : 0;
854 12 : break;
855 : case '>':
856 62 : v = l > r ? 1 : 0;
857 62 : break;
858 : case '&':
859 3 : v = l == 1 && r == 1 ? 1 : 0;
860 3 : break;
861 : case '|':
862 3 : v = l == 1 || r == 1 ? 1 : 0;
863 3 : break;
864 : case '=':
865 21 : v = l == r;
866 21 : break;
867 : case u'≠':
868 12 : v = l != r;
869 12 : break;
870 : case u'≤':
871 12 : v = l <= r ? 1 : 0;
872 12 : break;
873 : case u'≥':
874 12 : v = l >= r ? 1 : 0;
875 12 : break;
876 : case '^':
877 12 : v = pow(l, r);
878 12 : break;
879 : default:
880 0 : printf("unknown binary op (%c) encountered\n", n->op);
881 : }
882 16598 : break;
883 : case N_UNKNOWN:
884 : default:
885 0 : printf("unknown node encountered\n");
886 : // TODO: error
887 0 : break;
888 : }
889 :
890 2058940 : return v;
891 : }
892 :
893 : int
894 17 : sd_sim_get_stepcount(SDSim *sim)
895 : {
896 17 : if (!sim)
897 1 : return -1;
898 16 : return sim->nsaves;
899 : }
900 :
901 : int
902 15 : sd_sim_get_varcount(SDSim *sim)
903 : {
904 15 : if (!sim)
905 1 : return -1;
906 14 : return sim->nvars;
907 : }
908 :
909 : int
910 16 : sd_sim_get_varnames(SDSim *sim, const char **result, size_t max)
911 : {
912 16 : if (!sim || !result)
913 1 : return -1;
914 15 : else if (max == 0)
915 1 : return 0;
916 :
917 14 : return module_get_varnames(sim->module, result, max);
918 : }
919 :
920 : int
921 16 : module_get_varnames(AVar *module, const char **result, size_t max)
922 : {
923 16 : const char **start = result;
924 :
925 95 : for (size_t i = 0; i < module->avars.len && max > 0; i++) {
926 79 : AVar *av = module->avars.elems[i];
927 79 : if (av->model) {
928 2 : size_t n = module_get_varnames(av, result, max);
929 2 : result += n;
930 2 : max -= n;
931 77 : } else if (!av->src) {
932 : // only include non-ghosts in output
933 73 : *result = avar_qual_name(av);
934 73 : result++;
935 73 : max--;
936 : }
937 : }
938 16 : return result-start;
939 : }
940 :
941 : void
942 24 : module_clear_visited(AVar *module)
943 : {
944 172 : for (size_t i = 0; i < module->avars.len; i++) {
945 148 : AVar *av = module->avars.elems[i];
946 148 : if (av->model) {
947 4 : module_clear_visited(av);
948 : } else {
949 144 : av->visited = false;
950 144 : av->visiting = false;
951 : }
952 : }
953 24 : module->visited = false;
954 24 : module->visiting = false;
955 24 : }
956 :
957 : int
958 77 : sd_sim_get_series(SDSim *s, const char *name, double *results, size_t len)
959 : {
960 : int off;
961 : size_t i;
962 :
963 77 : if (!s || !name || !results)
964 1 : return -1;
965 :
966 76 : if (strcmp(name, "time") == 0) {
967 15 : off = 0;
968 : } else {
969 61 : AVar *av = resolve(s->module, name);
970 61 : if (!av)
971 1 : return -1;
972 60 : off = av->offset;
973 : }
974 :
975 32750 : for (i = 0; i <= s->nsaves && i < len; i++)
976 32675 : results[i] = s->slab[i*s->nvars + off];
977 :
978 75 : return i;
979 : }
980 :
981 : double
982 50 : rt_pulse(SDSim *s, Node *n, double dt, double time, size_t len, double *args)
983 : {
984 : double magnitude, first_pulse, next_pulse, interval;
985 :
986 50 : magnitude = args[0];
987 50 : first_pulse = args[1];
988 50 : if (len > 2)
989 50 : interval = args[2];
990 : else
991 0 : interval = 0;
992 :
993 50 : if (time < first_pulse)
994 16 : return 0;
995 :
996 34 : next_pulse = first_pulse;
997 :
998 68 : while (time >= next_pulse) {
999 34 : if (time < next_pulse + dt)
1000 2 : return magnitude/dt;
1001 32 : else if (interval <= 0)
1002 32 : break;
1003 : else
1004 0 : next_pulse += interval;
1005 : }
1006 32 : return 0;
1007 : }
1008 :
1009 : double
1010 12 : rt_min(SDSim *s, Node *n, double dt, double time, size_t len, double *args)
1011 : {
1012 : double a, b;
1013 :
1014 12 : if (len != 2)
1015 0 : return NAN;
1016 :
1017 12 : a = args[0];
1018 12 : b = args[1];
1019 :
1020 12 : return a < b ? a : b;
1021 : }
1022 :
1023 : double
1024 12 : rt_max(SDSim *s, Node *n, double dt, double time, size_t len, double *args)
1025 : {
1026 : double a, b;
1027 :
1028 12 : if (len != 2)
1029 0 : return NAN;
1030 :
1031 12 : a = args[0];
1032 12 : b = args[1];
1033 :
1034 12 : return a > b ? a : b;
1035 : }
|
