improved simulation stuff

This commit is contained in:
Jay Robson 2024-01-18 18:00:39 +11:00
parent e4f8391e39
commit 67773f14bd
19 changed files with 335 additions and 254 deletions

56
src/coolant/valve.hpp Normal file
View File

@ -0,0 +1,56 @@
#pragma once
namespace sim::coolant
{
template <class A>
class valve
{
A* a;
const double max;
const double pressure;
double state = 0;
double rate = 0;
public:
constexpr valve(A& a, double max, double pressure) : a(&a), max(max), pressure(pressure)
{
}
constexpr double get_state() const
{
return state;
}
constexpr void set_state(double v)
{
if(v > 1) v = 1;
if(v < 0) v = 0;
state = v;
}
constexpr void open(double v)
{
set_state(state + v);
}
constexpr void update(double secs)
{
rate = a->extract_steam(secs, state * max, pressure) / secs;
}
friend std::ostream& operator<<(std::ostream& o, const valve& v)
{
o << "Opened: " << (v.state * 100) << " %\n";
o << "Rate: " << v.rate << " g/s\n";
return o;
}
};
};

View File

@ -6,6 +6,7 @@
#include "reactor/coolant/heater.hpp"
#include "reactor/coolant/vessel.hpp"
#include "coolant/fluid_t.hpp"
#include "coolant/valve.hpp"
#include "display.hpp"
#include <cmath>
@ -14,6 +15,9 @@
#include <curses.h>
#include <sys/time.h>
const bool do_graphics = true;
const bool do_clock = true;
unsigned long get_now()
{
struct timeval tv;
@ -22,6 +26,11 @@ unsigned long get_now()
}
int main()
{
std::random_device rd;
std::mt19937 rand(rd());
if(do_graphics)
{
initscr();
cbreak();
@ -29,23 +38,27 @@ int main()
keypad(stdscr, TRUE);
nodelay(stdscr, TRUE);
curs_set(0);
}
sim::reactor::coolant::vessel vessel(200, 400, sim::coolant::WATER);
sim::reactor::reactor<5, 5> reactor = sim::reactor::builder<5, 5>(
sim::reactor::fuel::fuel_rod(2000, 4000),
sim::reactor::control::control_rod(10000, 1),
sim::reactor::control::control_rod(vessel, 10000, 1),
sim::reactor::coolant::pipe(vessel), {
"# ##",
"#FCF ",
" CHC ",
" FCF#",
"## #"
"#C#C#",
"CFCFC",
"#C#C#",
"CFCFC",
"#C#C#"
});
sim::coolant::valve<sim::reactor::coolant::vessel> valve(vessel, 1, 101000);
double secs = 0;
long clock = get_now();
double speed = 1;
int framerate = 100;
int steps_extra = 1;
for(;;)
{
@ -80,29 +93,23 @@ secs: ss << s << "s\n";
ss << "Speed: " << speed << "x\n\n";
}
int skip = 1;
while(speed / framerate / skip > 1)
for(int i = 0; i < steps_extra; i++)
{
skip *= 2;
}
for(int i = 0; i < skip; i++)
{
double dt = speed / framerate / skip;
reactor.update(dt);
double dt = speed / framerate / steps_extra;
reactor.update(rand, dt);
vessel.update();
vessel.extract_steam(dt, 0.01, 101000);
valve.update(dt);
secs += dt;
}
secs += speed / framerate;
ss << "Vessel\n" << vessel << "\n";
ss << "Steam Valve\n" << valve << "\n";
if(do_graphics)
{
erase();
display::draw_text(1, 0, ss.str().c_str());
}
const int X = 1, Y = 30;
const int W = 36, H = 11;
@ -110,7 +117,8 @@ secs: ss << s << "s\n";
for(int x = 0; x < reactor.width; x++)
for(int y = 0; y < reactor.height; y++)
{
sim::reactor::rod* r = reactor.rods[x][y];
int id = y * reactor.width + x;
sim::reactor::rod* r = reactor.rods[id];
if(!r->should_display())
{
@ -123,23 +131,30 @@ secs: ss << s << "s\n";
int px = X + (H - 1) * y;
int py = Y + (W - 1) * x;
if(do_graphics)
{
display::draw_text(px + 1, py + 2, ss.str().c_str());
display::draw_box(px, py, H, W);
}
if(r->should_select() && x == reactor.cursor_x && y == reactor.cursor_y)
if(do_graphics && r->should_select() && id == reactor.cursor)
{
display::draw_text(px + 1, py + W - 5, "[ ]");
}
if(r->is_selected())
if(do_graphics && r->is_selected())
{
display::draw_text(px + 1, py + W - 4, "#");
}
}
refresh();
int c = 0;
int c = getch();
if(do_graphics)
{
refresh();
c = getch();
}
switch(c)
{
@ -150,10 +165,10 @@ secs: ss << s << "s\n";
reactor.move_cursor(1);
break;
case KEY_UP:
reactor.update_selected(0.01);
reactor.update_selected(1);
break;
case KEY_DOWN:
reactor.update_selected(-0.01);
reactor.update_selected(-1);
break;
case ' ':
reactor.toggle_selected();
@ -164,8 +179,16 @@ secs: ss << s << "s\n";
case 'g':
speed /= 10;
break;
case 'r':
valve.open(0.001);
break;
case 'f':
valve.open(-0.001);
break;
}
if(do_clock)
{
long now = get_now();
while(clock + 1e6 / framerate > now)
@ -176,6 +199,7 @@ secs: ss << s << "s\n";
clock += 1e6 / framerate;
}
}
return 0;
}

View File

@ -5,10 +5,9 @@
using namespace sim::reactor::control;
control_rod::control_rod(double limit, double max)
control_rod::control_rod(coolant::vessel& v, double limit, double max) : coolant::pipe(v), limit(limit), max(max)
{
this->limit = limit;
this->max = max;
}
void control_rod::display(std::ostream& o) const
@ -17,11 +16,6 @@ void control_rod::display(std::ostream& o) const
o << "Use: " << absorbed << " / " << limit << " mol\n";
};
double control_rod::get_k(val_t type) const
{
return 0.5;
}
void control_rod::set_reactivity(double a)
{
inserted = 1 - a;
@ -31,13 +25,16 @@ void control_rod::update(double secs)
{
update_rod(secs);
double m = 1 - std::pow(0.5, (1 - absorbed / limit) * inserted * max);
double k = (1 - absorbed / limit) * inserted * max;
double m = 1 - std::pow(0.5, secs * -std::log2(1 - k));
double r_fast = vals[val_t::N_FAST] * m;
double r_slow = vals[val_t::N_SLOW] * m;
vals[val_t::N_FAST] -= r_fast;
vals[val_t::N_SLOW] -= r_slow;
absorbed += r_fast + r_slow;
update_pipe(secs);
}
void control_rod::update_selected(double a)

View File

@ -1,26 +1,26 @@
#pragma once
#include "../rod.hpp"
#include "../coolant/pipe.hpp"
namespace sim::reactor::control
{
class control_rod : public sim::reactor::rod
class control_rod : public sim::reactor::coolant::pipe
{
const double limit;
const double max;
double inserted = 1;
double absorbed = 0;
double limit;
double max;
virtual double get_k(sim::reactor::rod::val_t type) const;
virtual void display(std::ostream& o) const;
virtual const char* get_name() const { return "Control Rod"; }
public:
control_rod(double limit, double max);
control_rod(coolant::vessel& v, double limit, double max);
virtual void update(double secs);
void set_reactivity(double a);

View File

@ -26,9 +26,8 @@ void graphite_rod::update(double secs)
{
update_rod(secs);
double v = vals[val_t::N_FAST];
vals_in[val_t::N_FAST] -= v;
vals_in[val_t::N_SLOW] += v;
vals[val_t::N_SLOW] = vals[val_t::N_FAST];
vals[val_t::N_FAST] = 0;
}
void graphite_rod::update_selected(double a)

View File

@ -10,10 +10,10 @@ class graphite_rod : public sim::reactor::rod
{
double inserted = 0;
virtual double get_k(sim::reactor::rod::val_t type) const;
virtual void display(std::ostream& o) const;
virtual const char* get_name() const { return "Graphite Rod"; }
virtual double get_k(val_t type) const;
public:

View File

@ -7,7 +7,7 @@ void heater::update(double secs)
{
update_rod(secs);
vals_in[val_t::HEAT] += rate * secs;
vals[val_t::HEAT] += rate * secs;
}
void heater::update_selected(double a)
@ -15,11 +15,6 @@ void heater::update_selected(double a)
rate += a;
}
double heater::get_k(sim::reactor::rod::val_t type) const
{
return 0.5;
}
void heater::display(std::ostream& o) const
{
o << "Rate: " << rate << " C/s\n";

View File

@ -10,10 +10,10 @@ class heater : public sim::reactor::rod
{
double rate = 0;
virtual double get_k(sim::reactor::rod::val_t type) const;
virtual void display(std::ostream& o) const;
virtual const char* get_name() const { return "Heater"; }
virtual double get_k(val_t type) const { return 0.5; }
public:

View File

@ -16,16 +16,14 @@ double pipe::get_k(val_t type) const
void pipe::update(double secs)
{
double v;
update_rod(secs);
v = vessel->add_heat(vals[val_t::HEAT]);
steam = vals[val_t::HEAT] - v;
vals[val_t::HEAT] = v;
v = vals[val_t::N_FAST];
vals_in[val_t::N_FAST] -= v;
vals_in[val_t::N_SLOW] += v;
update_pipe(secs);
}
void pipe::update_pipe(double secs)
{
vals[val_t::HEAT] = vessel->add_heat(vals[val_t::HEAT]);
vals[val_t::N_SLOW] += vals[val_t::N_FAST];
vals[val_t::N_FAST] = 0;
}

View File

@ -9,20 +9,22 @@ namespace sim::reactor::coolant
class pipe : public sim::reactor::rod
{
protected:
coolant::vessel* vessel;
double steam;
virtual double get_k(sim::reactor::rod::val_t type) const;
virtual const char* get_name() const { return "Coolant"; }
void update_pipe(double secs);
public:
pipe(coolant::vessel& v);
virtual bool should_display() const { return true; }
virtual void update(double secs);
};
}

View File

@ -21,6 +21,7 @@ void vessel::update()
double s = steam + n;
double l = fluid.l_to_g(level) - n;
double v = fluid.l_to_g(volume);
if(l < 0)
{
@ -28,12 +29,18 @@ void vessel::update()
l = 0;
}
if(fluid.g_to_l(s) > volume)
if(s > v)
{
s = fluid.l_to_g(volume);
s = v;
l = 0;
}
if(l > v)
{
l = v;
s = 0;
}
double diff = s - steam;
steam = s;
@ -56,33 +63,22 @@ double vessel::extract_steam(double dt, double a, double p2)
{
// calculate the mass moved
double p1 = get_pressure();
double p = p1 - p2;
double m = 1;
double p = (p1 - p2) * 0.001; // mPa or g/m/s^2
if(p == 0)
{
return 0;
}
if(p < 0)
{
m = -1;
p = -p;
double V = (volume - level) * 0.001; // m^3
double mass = std::min(dt * a * p / std::sqrt( V * std::abs(p) / steam ), steam);
if(std::isnan(mass))
{
return 0;
}
double V = (volume - level) * 0.001;
double v = std::sqrt( V * p / steam );
double mass = m * dt * a * p / v;
if(mass > steam)
{
mass = steam;
}
steam -= mass;
return mass;
}

View File

@ -12,26 +12,22 @@ void fuel_rod::display(std::ostream& o) const
{
o << "Fuel: " << s.get_fuel() << " / " << s.get_mass() << " mol\n";
o << "Efficiency: " << (s.get_efficiency() * 100) << " %\n";
o << "Energy: +" << s.get_energy() << " C\n";
o << "Output: " << s.get_energy() << " W/s\n";
o << "Iodine: " << s.get_i_135() << " mol\n";
o << "Xenon: " << s.get_xe_135() << " mol\n";
}
double fuel_rod::get_k(val_t type) const
{
return 0.5;
}
void fuel_rod::update(double secs)
{
update_rod(secs);
s.clear_energy();
s.clear_fast_neutrons();
s.add_slow_neutrons(vals[val_t::N_SLOW]);
vals_in[val_t::HEAT] += s.extract_energy();
vals_in[val_t::N_FAST] += s.extract_fast_neutrons();
vals_in[val_t::N_SLOW] -= vals[val_t::N_SLOW];
s.update(secs);
vals[val_t::HEAT] += s.get_energy() * secs;
vals[val_t::N_FAST] += s.get_fast_neutrons();
vals[val_t::N_SLOW] = 0;
}

View File

@ -11,7 +11,7 @@ class fuel_rod : public sim::reactor::rod
{
sample s;
virtual double get_k(sim::reactor::rod::val_t type) const;
virtual double get_k(val_t type) const { return 0.5; }
virtual void display(std::ostream& o) const;
virtual const char* get_name() const { return "Fuel"; }

View File

@ -2,6 +2,8 @@
#include "sample.hpp"
#include "half_life.hpp"
#include <iostream>
using namespace sim::reactor::fuel;
static const double NEUTRON_BG = 1e-10;
@ -12,44 +14,11 @@ sample::sample(double fuel, double mass)
this->mass = mass;
}
void sample::absorb_slow_neutrons()
{
// absorb neutrons
double volume = get_volume();
double neutrons = slow_neutrons;
double neutrons_fuel = neutrons * (fuel / volume);
double neutrons_iodine = neutrons * (i_135 / volume);
double neutrons_xenon = neutrons * ((xe_135 * Xe_135_M) / volume);
double neutrons_total = neutrons;
slow_neutrons = 0;
xe_135 -= neutrons_xenon;
i_135 -= neutrons_iodine;
fuel -= neutrons_fuel;
// do the poison
te_135 += neutrons_fuel * (1.0 / 50.0);
xe_135 -= neutrons_xenon;
i_135 -= neutrons_iodine;
// deal with these edge cases
if(xe_135 < 0) { slow_neutrons -= xe_135; xe_135 = 0; }
if(i_135 < 0) { slow_neutrons -= i_135; i_135 = 0; }
if(fuel < 0) { slow_neutrons -= fuel; fuel = 0; }
efficiency = neutrons_fuel / neutrons_total;
// simulate fuel use
energy += neutrons_fuel;
waste.add_fissile(neutrons_fuel * 6);
}
void sample::update(double secs)
{
double m;
// sim::reactor::fuelulate waste and extract products
// decay waste and extract products
waste.update(secs);
fast_neutrons += waste.extract_neutrons();
energy += waste.extract_energy();
@ -69,27 +38,29 @@ void sample::update(double secs)
te_135 *= m;
// absorb neutrons
slow_neutrons += NEUTRON_BG * secs;
absorb_slow_neutrons();
}
double volume = get_volume();
double neutrons = slow_neutrons + NEUTRON_BG * secs;
double neutrons_fuel = neutrons * (fuel / volume);
double neutrons_iodine = neutrons * (i_135 / volume);
double neutrons_xenon = neutrons * ((xe_135 * Xe_135_M) / volume);
double neutrons_total = neutrons;
slow_neutrons = 0;
double sample::extract_energy()
{
double v = energy;
energy = 0;
return v;
}
te_135 += neutrons_fuel * (1.0 / 50.0);
xe_135 -= neutrons_xenon;
i_135 -= neutrons_iodine;
fuel -= neutrons_fuel;
double sample::extract_fast_neutrons()
{
double v = fast_neutrons;
fast_neutrons = 0;
return v;
}
// deal with these edge cases
if(xe_135 < 0) { slow_neutrons += xe_135; xe_135 = 0; }
if(i_135 < 0) { slow_neutrons += i_135; i_135 = 0; }
if(fuel < 0) { slow_neutrons += fuel; fuel = 0; }
void sample::add_slow_neutrons(double a)
{
slow_neutrons += a;
efficiency = neutrons_fuel / neutrons_total;
// simulate fuel use
energy += neutrons_fuel / secs;
waste.add_fissile(neutrons_fuel * 6);
}
void sample::display(std::ostream& o) const

View File

@ -14,38 +14,41 @@ class sample
sim::reactor::fuel::waste waste;
// mol
double fuel = 0;
double i_135 = 0;
double xe_135 = 0;
double te_135 = 0;
double mass = 0;
double energy = 0;
double fast_neutrons = 0;
double slow_neutrons = 0;
double energy = 0; // W/s
double fast_neutrons = 0; // mol/s
double slow_neutrons = 0; // mol/s
double efficiency = 0;
void display(std::ostream& o) const;
void absorb_slow_neutrons();
void absorb_slow_neutrons(double secs);
public:
sample(double fuel, double mass);
void update(double secs);
double extract_energy();
double extract_fast_neutrons();
void add_slow_neutrons(double a);
constexpr double get_fuel() const { return fuel; }
constexpr double get_mass() const { return mass; }
constexpr double get_energy() const { return energy; }
constexpr double get_fast_neutrons() const { return fast_neutrons; }
constexpr double get_volume() const { return mass + xe_135 * Xe_135_M; }
constexpr double get_efficiency() const { return efficiency; }
constexpr double get_te_135() const { return te_135; }
constexpr double get_i_135() const { return i_135; }
constexpr double get_xe_135() const { return xe_135; }
constexpr void clear_energy() { energy = 0; }
constexpr void clear_fast_neutrons() { fast_neutrons = 0; }
constexpr void add_slow_neutrons(double a) { slow_neutrons += a; }
friend std::ostream& operator<<(std::ostream& o, const sample& s)
{
s.display(o);

View File

@ -6,10 +6,10 @@ using namespace sim::reactor::fuel;
void waste::update(double secs)
{
double hl = 1;
double next[waste::N - 1] = {0};
double hl = 1;
for(int i = 0; i < waste::N - 1; i++)
for(int i = 0; i < waste::N - 1; hl *= 2, i++)
{
double m = 1 - half_life::get(secs, hl);
double h = high[i] * m;
@ -21,7 +21,6 @@ void waste::update(double secs)
neutrons += h;
energy += h + l;
hl *= 2;
}
for(int i = 0; i < waste::N - 1; i++)

View File

@ -4,6 +4,9 @@
#include "rod.hpp"
#include <array>
#include <random>
#include <algorithm>
#include <iostream>
namespace sim::reactor
{
@ -11,99 +14,113 @@ namespace sim::reactor
template <int W, int H>
struct reactor
{
const static int width = W;
const static int height = H;
constexpr const static int width = W;
constexpr const static int height = H;
constexpr const static int size = W*H;
std::array<std::array<rod*, H>, W> rods;
rod* rods[size];
int cursor_x = 0;
int cursor_y = 0;
int cursor = 0;
reactor(std::array<rod*, W * H> rods)
reactor(std::array<rod*, size> rods)
{
for(int y = 0; y < H; y++)
for(int x = 0; x < W; x++)
for(int i = 0; i < size; i++)
{
this->rods[x][y] = rods[y * W + x];
this->rods[i] = rods[i];
}
}
void update(double secs)
void update(std::mt19937& rand, double secs)
{
// do interactions
for(int x = 1; x < W; x++)
int rods_lookup[size];
for(int i = 0; i < size; i++)
{
rods[x][0]->interact(rods[x - 1][0], secs);
rods_lookup[i] = i;
}
for(int y = 1; y < H; y++)
for(int i = 0; i < size; i++)
{
rods[0][y]->interact(rods[0][y - 1], secs);
rods[i]->update(secs);
}
for(int y = 1; y < H; y++)
for(int x = 1; x < W; x++)
{
rod* r = rods[x][y];
r->interact(rods[x - 1][y], secs);
r->interact(rods[x][y - 1], secs);
update_interactions(rand, rods_lookup, secs / 2);
}
// do updates
for(int y = 0; y < H; y++)
for(int x = 0; x < W; x++)
void update_selected(int v)
{
rods[x][y]->update(secs);
for(int i = 0; i < size; i++)
{
rod* r = rods[i];
if(r->is_selected())
{
r->update_rod_selected(v);
}
}
}
void update_selected(double v)
int move_cursor(int d)
{
for(int y = 0; y < H; y++)
for(int x = 0; x < W; x++)
if(rods[x][y]->is_selected())
for(int i = 0; i < size; i++)
{
rods[x][y]->update_selected(v);
cursor = (cursor + d) % size;
if(cursor < 0)
{
cursor += size;
}
if(rods[cursor]->should_select())
{
return cursor;
}
}
void move_cursor(int d)
{
for(;;)
{
cursor_x += d;
while(cursor_x >= W)
{
cursor_x -= W;
cursor_y += 1;
}
while(cursor_x < 0)
{
cursor_x += W;
cursor_y -= 1;
}
cursor_y %= H;
if(cursor_y < 0)
{
cursor_y += H;
}
if(rods[cursor_x][cursor_y]->should_select())
{
return;
}
}
return 0;
}
void toggle_selected()
{
if(rods[cursor_x][cursor_y]->should_select())
if(rods[cursor]->should_select())
{
rods[cursor_x][cursor_y]->toggle_selected();
rods[cursor]->toggle_selected();
}
}
private:
void update_tile(std::mt19937& rand, double secs, int i, int x, int y)
{
std::array<int, 2> nb_lookup[4] = {{-1, 0}, {1, 0}, {0, -1}, {0, 1}};
std::shuffle(nb_lookup, &nb_lookup[3], rand);
for(int j = 0; j < 4; j++)
{
int xp = x + nb_lookup[j][0];
int yp = y + nb_lookup[j][1];
if(xp >= 0 && yp >= 0 && xp < width && yp < height)
{
rods[i]->interact(rods[yp * width + xp], secs / 2);
}
}
}
void update_interactions(std::mt19937& rand, int* rods_lookup, double secs)
{
std::shuffle(rods_lookup, &rods_lookup[size - 1], rand);
for(int id = 0; id < size; id++)
{
int i = rods_lookup[id];
int x = i % width;
int y = i / width;
for(int j = 0; j < 4; j++)
{
update_tile(rand, secs, i, x, y);
}
}
}
};

View File

@ -12,7 +12,7 @@ double rod::get(val_t type) const
void rod::add(val_t type, double v)
{
vals_in[type] += v;
vals[type] += v;
}
double rod::extract(val_t type, double s, double k, double o)
@ -20,14 +20,16 @@ double rod::extract(val_t type, double s, double k, double o)
k *= get_k(type);
double m = 1;
k = 1 - k * get_k(type);
if(k < 1)
if(k > 0)
{
m = 1 - std::pow(0.5, s * -std::log2(1 - k * get_k(type)));
m = 1 - std::pow(0.5, s * -std::log2(k));
}
double v = m * 0.5 * (get(type) - o);
vals_in[type] -= v;
vals[type] -= v;
return v;
}
@ -40,18 +42,35 @@ void rod::interact(rod* o, double secs)
}
}
void rod::update_rod(double secs)
double rod::get_speed() const
{
for(int i = 0; i < rod::VAL_N; i++)
{
val_t v = (val_t)i;
vals[v] += vals_in[v];
vals_in[v] = 0;
int m = motion < 0 ? -1 : 1;
return motion == 0 ? 0 : (std::pow(10, std::abs(motion)) * 1e-10 * m);
}
void rod::update_rod(double secs)
{
// decay the free neutrons
double m = std::pow(0.5, secs / 879.4);
vals[val_t::N_FAST] *= m;
vals[val_t::N_SLOW] *= m;
if(motion != 0 && !is_selected())
{
motion = 0;
}
if(motion != 0)
{
update_selected(get_speed() * secs);
}
}
void rod::update_rod_selected(int m)
{
motion += m;
if(motion > 10) motion = 10;
if(motion < -10) motion = -10;
}

View File

@ -27,17 +27,24 @@ public:
virtual bool should_display() const { return false; }
virtual bool should_select() const { return false; }
virtual void update_selected(double a) { }
void update_rod_selected(int m);
constexpr void toggle_selected() { selected = !selected; }
constexpr bool is_selected() { return selected; }
constexpr bool is_selected() const { return selected; }
friend std::ostream& operator<<(std::ostream& o, const rod& r)
{
if(!r.should_display()) return o;
o << r.get_name() << "\n";
if(r.is_selected())
{
o << "Speed: " << r.get_speed() << "\n";
}
r.display(o);
o << "Heat: " << r.get(val_t::HEAT) << " C\n";
o << "Fast: " << r.get(val_t::N_FAST) << " mol\n";
o << "Slow: " << r.get(val_t::N_SLOW) << " mol\n";
@ -47,15 +54,17 @@ public:
protected:
double vals_in[VAL_N] = {0};
double vals[VAL_N] = {0};
bool selected = false;
int motion = 0;
virtual void display(std::ostream& o) const { };
virtual double get_k(val_t type) const { return 0; }
virtual const char* get_name() const { return "Empty"; }
virtual void update_selected(double a) { }
void update_rod(double secs);
double get_speed() const;
};
}