fdtd-coremem/crates/applications/archive/toroid25d.rs

use coremem::{Driver, mat, meas, MaterialSim as _, SimState};
use coremem::geom::{CylinderZ, Index, Meters, Vec2, Vec3};
use coremem::real::R64 as Real;
use coremem::stim::{Stimulus, Sinusoid3};
use log::trace;

fn main() {
    coremem::init_logging();
    for p in 0..20 {
        let ferro_depth = 10e-6 * (20-p) as f32;
        let feat_size = 10e-6f32; // feature size
        let from_m = |m| (m/feat_size) as u32;
        let m_to_um = |px| (px * 1e6) as u32;
        let to_m = |px| px as f32 * feat_size;
        let width = 2500e-6;
        let depth = 200e-6;
        let conductor_inner_rad = 0e-6;
        let conductor_outer_rad = 19e-6;
        let ferro_inner_rad = 100e-6;
        let ferro_outer_rad = 200e-6;
        //let ferro_depth = 10e-6;
        let buffer = 250e-6;
        let peak_current = 2e3;
        let current_duration = 1e-9; // half-wavelength of the sine wave
        let conductivity = 1.0e5;
        let half_width = width * 0.5;
        let half_depth = depth * 0.5;
        let steps_per_frame = 160;

        let duration = 1e-9;

        let width_px = from_m(width);
        let depth_px = from_m(depth);
        let size_px = Index((width_px, width_px, depth_px).into());
        let mut driver: Driver<SimState<Real>> = Driver::new(size_px, feat_size);
        let base = "toroid25d-9";
        let _ = std::fs::create_dir(base);
        let prefix = format!("{}/{}-flt{}-{}-feat{}um-{}mA-{}ps--radii{}um-{}um-{}um-{}um",
            base,
            base,
            std::mem::size_of::<Real>() * 8,
            *size_px,
            m_to_um(feat_size),
            (peak_current * 1e3) as i64,
            (current_duration * 1e12) as i64,
            m_to_um(conductor_outer_rad),
            m_to_um(ferro_inner_rad),
            m_to_um(ferro_outer_rad),
            m_to_um(ferro_depth),
        );
        let _ = std::fs::create_dir(&prefix);
        //driver.add_y4m_renderer(&*format!("{}.y4m", prefix), steps_per_frame);
        //driver.add_plotly_renderer(&*format!("{}/frame-", prefix), steps_per_frame);
        driver.add_serializer_renderer(&*format!("{}/frame-", prefix), steps_per_frame, None);
        let conductor_region = CylinderZ::new(
            Vec2::new(half_width, half_width),
            conductor_outer_rad);
        // driver.add_term_renderer(steps_per_frame);
        driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current)));
        //driver.add_measurement(meas::Current(conductor_region.clone()));
        driver.add_measurement(meas::MagnetizationAt(
            Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
        ));
        driver.add_measurement(meas::MagneticFluxAt(
            Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
        ));
        driver.add_measurement(meas::MagneticStrengthAt(
            Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
        ));
        driver.add_measurement(meas::MagnetizationAt(
            Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
        ));
        driver.add_measurement(meas::MagneticFluxAt(
            Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
        ));
        driver.add_measurement(meas::MagneticStrengthAt(
            Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
        ));
        driver.add_measurement(meas::MagnetizationAt(
            Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
        ));
        driver.add_measurement(meas::MagneticFluxAt(
            Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
        ));
        driver.add_measurement(meas::MagneticStrengthAt(
            Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
        ));

        let center = Vec2::new(half_width, half_width);

        for z_px in 0..depth_px {
            for y_px in 0..width_px {
                for x_px in 0..width_px { 
                    let loc = Index((x_px, y_px, z_px).into());
                    let d = Vec2::new(to_m(x_px), to_m(y_px)) - center;
                    let r = d.mag();
                    if (conductor_inner_rad..conductor_outer_rad).contains(&r) {
                        driver.state.put_material(loc, mat::IsomorphicConductor::new(conductivity));
                    } else if (ferro_inner_rad..ferro_outer_rad).contains(&r) {
                        let half_depth_px = from_m(half_depth);
                        let ferro_depth_px = from_m(ferro_depth);
                        if (half_depth_px-ferro_depth_px/2 .. half_depth_px+(ferro_depth_px+1)/2).contains(&z_px) {
                            trace!("placing ferro at {:?}", loc);
                            driver.state.put_material(loc, mat::db::ferroxcube_3r1());
                        }
                    }
                }
            }
        }
        let boundary_xy = half_width - ferro_outer_rad - buffer;
        println!("boundary: {}um", m_to_um(boundary_xy));
        let boundary = Index((from_m(boundary_xy), from_m(boundary_xy), 0).into());
        driver.add_pml_boundary(boundary);

        driver.add_stimulus(Stimulus::new(
            conductor_region.clone(),
            Sinusoid3::from_wavelength(Vec3::new(0.0, 0.0, peak_current), current_duration * 2.0
    )));

        while driver.dyn_state().time() < duration {
            driver.step();
        }
    }
}