WIP Stimulus system

examples/toroid25d.rs

@@ -1,5 +1,6 @@
 use coremem::{consts, Driver, Flt, mat, meas};
-use coremem::geom::{Coord, CylinderZ, Vec2};
+use coremem::geom::{Coord, CylinderZ, Vec2, Vec3};
+use coremem::stim::{Stimulus, Sinusoid};
 
 fn main() {
     coremem::init_logging();
@@ -35,12 +36,12 @@ fn main() {
         m_to_um(ferro_inner_rad),
         m_to_um(ferro_outer_rad),
     ));
+    let conductor_region = CylinderZ::new(
+        Vec2::new(half_width, half_width),
+        conductor_outer_rad);
     // driver.add_term_renderer();
     driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current)));
-    driver.add_measurement(meas::Current(CylinderZ::new(
+    driver.add_measurement(meas::Current(conductor_region.clone()));
-        Vec2::new(half_width, half_width),
-        conductor_outer_rad)
-    ));
     driver.add_measurement(meas::Magnetization(
         (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()
     ));
@@ -92,28 +93,11 @@ fn main() {
     let boundary = Coord::new(from_m(boundary_xy), from_m(boundary_xy), 20);
     driver.add_upml_boundary(boundary);
 
+    driver.add_stimulus(Stimulus::new(
+        conductor_region.clone(),
+        Sinusoid::new(Vec3::new(0.0, 0.0, peak_current * 1e-18), 1e9)));
+
     loop {
-        // let drive_current = peak_current * match driver.state.step_no() {
-        //     0..=1000 => 1.0,
-        //     3000..=4000 => -1.0,
-        //     _ => 0.0,
-        // };
-        let drive_current = peak_current * 1e-18;
-        // J = \sigma*E = [Am^-2]
-        // I = \sigma*E*Area
-        // E = I / \sigma / Area
-        //let e = v/(2.0*feat_size);
-        let area = consts::PI*(conductor_outer_rad*conductor_outer_rad - conductor_inner_rad*conductor_inner_rad);
-        let e = drive_current/conductivity/area;
-        for y_px in from_m(half_width-conductor_outer_rad)..from_m(half_width+conductor_outer_rad) {
-            for x_px in from_m(half_width-conductor_outer_rad)..from_m(half_width+conductor_outer_rad) { 
-                let loc = Coord::new(x_px, y_px, 0);
-                let d = Vec2::new(to_m(x_px), to_m(y_px)) - center;
-                if (conductor_inner_rad..conductor_outer_rad).contains(&d.mag()) {
-                    driver.state.impulse_ez(loc, e);
-                }
-            }
-        }
         driver.step();
     }
 }

src/driver.rs

@@ -4,6 +4,7 @@ use crate::geom::{Coord, Vec3};
 use crate::meas::{self, AbstractMeasurement};
 use crate::render::{self, MultiRenderer, Renderer};
 use crate::sim::{GenericSim as _, SimState};
+use crate::stim::AbstractStimulus;
 
 use log::{info, trace};
 use std::path::PathBuf;
@@ -14,25 +15,32 @@ pub struct Driver {
     renderer: MultiRenderer,
     steps_per_frame: u64,
     time_spent_stepping: Duration,
+    time_spent_on_stimuli: Duration,
     time_spent_rendering: Duration,
     measurements: Vec<Box<dyn AbstractMeasurement>>,
+    stimuli: Vec<Box<dyn AbstractStimulus>>,
     start_time: SystemTime,
 }
 
 impl Driver {
-    // TODO: allow depth
     pub fn new(size: Coord, feature_size: Flt) -> Self {
         Driver {
             state: SimState::new(size, feature_size),
             renderer: Default::default(),
             steps_per_frame: 1,
             time_spent_stepping: Default::default(),
+            time_spent_on_stimuli: Default::default(),
             time_spent_rendering: Default::default(),
             measurements: vec![Box::new(meas::Time), Box::new(meas::Meta), Box::new(meas::Energy)],
+            stimuli: vec![],
             start_time: SystemTime::now(),
         }
     }
 
+    pub fn add_stimulus<S: AbstractStimulus + 'static>(&mut self, s: S) {
+        self.stimuli.push(Box::new(s))
+    }
+
     pub fn add_measurement<M: AbstractMeasurement + 'static>(&mut self, m: M) {
         self.measurements.push(Box::new(m));
     }
@@ -126,6 +134,15 @@ impl Driver {
             self.time_spent_rendering += start_time.elapsed().unwrap();
             trace!("render end");
         }
+        {
+            trace!("stimuli begin");
+            let start_time = SystemTime::now();
+            for stim in &mut *self.stimuli {
+                stim.apply(&mut self.state);
+            }
+            self.time_spent_on_stimuli += start_time.elapsed().unwrap();
+        }
+
         trace!("step begin");
         let start_time = SystemTime::now();
         self.state.step();
@@ -133,19 +150,21 @@ impl Driver {
         trace!("step end");
         let step = self.state.step_no();
         if step % (10*self.steps_per_frame) == 0 {
-            let driver_time = self.time_spent_stepping.as_secs_f64() as Flt;
+            let step_time = self.time_spent_stepping.as_secs_f64();
-            let render_time = self.time_spent_rendering.as_secs_f64() as Flt;
+            let stim_time = self.time_spent_on_stimuli.as_secs_f64();
-            let overall_time = self.start_time.elapsed().unwrap().as_secs_f64() as Flt;
+            let render_time = self.time_spent_rendering.as_secs_f64();
-            let fps = (self.state.step_no() as Flt) / overall_time;
+            let overall_time = self.start_time.elapsed().unwrap().as_secs_f64();
-            let sim_time = self.state.time();
+            let fps = (self.state.step_no() as f64) / overall_time;
+            let sim_time = self.state.time() as f64;
             info!(
-                "t={:.2e} frame {:06} fps: {:6.2} (sim: {:.1}s, render: {:.1}s, other: {:.1}s)",
+                "t={:.2e} frame {:06} fps: {:6.2} (sim: {:.1}s, stim: {:.1}s, render: {:.1}s, other: {:.1}s)",
                 sim_time,
                 step,
                 fps,
-                driver_time,
+                step_time,
+                stim_time,
                 render_time,
-                overall_time - driver_time - render_time
+                overall_time - step_time - stim_time - render_time
             );
         }
     }

src/geom/mod.rs

@@ -183,6 +183,7 @@ pub trait Region {
     fn contains(&self, p: Vec3) -> bool;
 }
 
+#[derive(Copy, Clone)]
 pub struct CylinderZ {
     center: Vec2,
     radius: Real,

src/lib.rs

@@ -13,6 +13,7 @@ pub mod mat;
 pub mod meas;
 pub mod render;
 pub mod sim;
+pub mod stim;
 
 pub use driver::*;
 pub use mat::*;
@@ -65,6 +66,7 @@ pub mod consts {
     // Vacuum Permittivity
     pub const EPS0: Flt = 8.854187812813e-12; // F⋅m−1
     pub const PI: Flt = std::f64::consts::PI as Flt;
+    pub const TWO_PI: Flt = 2.0 * std::f64::consts::PI as Flt;
     pub mod real {
         use super::*;
         pub(crate) fn C() -> Real {

src/stim.rs

@@ -0,0 +1,56 @@
+use crate::consts;
+use crate::flt::Flt;
+use crate::geom::{Region, Vec3};
+use crate::sim::GenericSim;
+
+pub trait AbstractStimulus {
+    // /// Returns a Region over which to apply some net impulse (E).
+    // /// i.e. the caller will divide E by the volume of the Region and then
+    // /// apply that to each cell in the region.
+    // fn e_for(&mut self, sim: &dyn GenericSim) -> (Box<Region>, Vec3);
+    fn apply(&mut self, sim: &mut dyn GenericSim);
+}
+
+pub trait TimeVarying {
+    /// Retrieve the E impulse to apply at the provided time (in seconds).
+    fn at(&mut self, t_sec: Flt) -> Vec3;
+}
+
+pub struct Stimulus<R, T> {
+    region: R,
+    time: T,
+}
+
+impl<R, T> Stimulus<R, T> {
+    pub fn new(region: R, time: T) -> Self {
+        Self {
+            region, time
+        }
+    }
+}
+
+impl<R: Region, T: TimeVarying> AbstractStimulus for Stimulus<R, T> {
+    fn apply(&mut self, sim: &mut dyn GenericSim) {
+        unimplemented!()
+    }
+}
+
+pub struct Sinusoid {
+    amp: Vec3,
+    omega: Flt,
+}
+
+impl Sinusoid {
+    pub fn new(amp: Vec3, freq: Flt) -> Self {
+        Self {
+            amp,
+            omega: freq * consts::TWO_PI,
+        }
+    }
+}
+
+impl TimeVarying for Sinusoid {
+    fn at(&mut self, t_sec: Flt) -> Vec3 {
+        self.amp * (t_sec * self.omega).sin()
+    }
+}