Rough working implementation of MPI parallelized synchronous reluctance machine optimization

apps/Switched Winding Synchrel/CMakeLists.txt

@@ -3,7 +3,10 @@ PROJECT(Switched\ Winding\ Synchrel)
 set(SOURCE_FILES
         main.cpp
         designspace.h designspace.cpp
-        globalobjective.h globalobjective.cpp)
+        globalobjective.h globalobjective.cpp
+        parameters.h parameters.cpp
+        model.h model.cpp
+        torquespeedcurve.h torquespeedcurve.cpp)
 
 add_executable(sws ${SOURCE_FILES})
 target_link_libraries(sws ${OERSTED_LIBRARIES})
\ No newline at end of file

apps/Switched Winding Synchrel/designspace.cpp

 #include "designspace.h"
 
+void limit_angular_thickness(CoordinateSpace &d, Particle &p) {
+    double_t sum_at{0.0};
+    for (size_t i = 0; i != 6; ++i) {
+        std::string key_at = "at" + std::to_string(i);
+        sum_at += p[key_at].Position;
+    }
+
+    if (sum_at > 0.93) {
+        for (size_t iter = 0; iter != 3; ++iter) {
+            double_t scale{0.93 / sum_at};
+            sum_at = 0;
+            for (size_t i = 0; i != 10; ++i) {
+                std::string key = "at" + std::to_string(i);
+
+                if (iter == 0) {
+                    p[key].Velocity = p[key].Position * (scale - 1.0);
+                }
+                p[key].Position = p[key].Position * scale;
+
+                if (p[key].Position < d[key].lower_bound()) {
+                    if (iter == 0) {
+                        p[key].Velocity = d[key].lower_bound() - p[key].Position;
+                    }
+                    p[key].Position = d[key].lower_bound();
+                }
+
+                sum_at += p[key].Position;
+            }
+        }
+    }
+};
+
+void limit_radial_thickness(CoordinateSpace &d, Particle &p) {
+    double_t sum_rt{0.0};
+    for (size_t i = 0; i != 7; ++i) {
+        std::string key = "rt" + std::to_string(i);
+        sum_rt += p[key].Position;
+    }
+
+    if (sum_rt > 1.0) {
+        for (size_t iter = 0; iter != 3; ++iter) {
+            double_t scale{1.0 / sum_rt};
+            sum_rt = 0;
+            for (size_t i = 0; i != 7; ++i) {
+                std::string key = "rt" + std::to_string(i);
+
+                if (iter == 0) {
+                    p[key].Velocity = p[key].Position * (scale - 1.0);
+                }
+                p[key].Position = p[key].Position * scale;
+
+                if (p[key].Position < d[key].lower_bound()) {
+                    if (iter == 0) {
+                        p[key].Velocity = d[key].lower_bound() - p[key].Position;
+                    }
+                    p[key].Position = d[key].lower_bound();
+                }
+
+                sum_rt += p[key].Position;
+            }
+        }
+    }
+};
+
 CoordinateSpace design_space() {
     CoordinateSpace d;
 
@@ -22,5 +86,8 @@ CoordinateSpace design_space() {
     d["bi"] = Coordinate{0.5,1.0,1.5};
     d["sd"] = Coordinate{2,6,14};
 
+    d.BoundaryProjections.push_back(limit_angular_thickness);
+    d.BoundaryProjections.push_back(limit_radial_thickness);
+
     return d;
-}
\ No newline at end of file
+}

apps/Switched Winding Synchrel/designspace.h

@@ -5,4 +5,8 @@
 
 CoordinateSpace design_space();
 
+void limit_angular_thicknes(CoordinateSpace &d, Particle &p);
+
+void limit_radial_thickenss(CoordinateSpace &d, Particle &p);
+
 #endif //OE_APP_SWITCHED_WINDING_SYNCHREL_DESIGNSPACE_H
\ No newline at end of file

apps/Switched Winding Synchrel/globalobjective.cpp

 #include "globalobjective.h"
+#include "model.h"
+#include "torquespeedcurve.h"
 
-ObjectiveMap evaluate(Particle p) {
-    // TODO
-    return ObjectiveMap{};
+//#include "parameters.h"
+//#include "create_sketch.h"
+//#include "create_mesh.h"
+//#include "create_physics.h"
+//#include "create_post_processor.h"
+
+ObjectiveMap evaluate(Particle particle) {
+    Model model{particle};
+
+    // Setup OpenMP
+    #ifdef _OPENMP
+        omp_set_num_threads(omp_get_max_threads());
+    #endif
+
+    size_t Na{5}, Nj{5};
+    Eigen::MatrixXd angle = Eigen::MatrixXd::Zero(Na,Nj);
+    Eigen::MatrixXd current = Eigen::MatrixXd::Zero(Na,Nj);
+    Eigen::MatrixXd torque = Eigen::MatrixXd::Zero(Na,Nj);
+    std::array<Eigen::MatrixXd,3> flux_cos;
+    std::array<Eigen::MatrixXd,3> flux_sin;
+    for (size_t i = 0; i != 3; ++i) {
+        flux_cos[i] = Eigen::MatrixXd::Zero(Na,Nj);
+        flux_sin[i] = Eigen::MatrixXd::Zero(Na,Nj);
+    }
+
+    std::vector<Oe::VectorXd> previous(Na * Nj);
+
+    for (size_t i = 0; i != Na; ++i) {
+        for (size_t j = 0; j != Nj; ++ j) {
+            angle(i,j) = i * M_PI_2 / (Na - 1);
+            current(i,j) = j * model.Params.Jmax / (Na - 1) * M_SQRT2 * model.Params.SlotFill;
+        }
+    }
+
+    // Torque-speed curve values
+    Eigen::ArrayXd SPEED = Eigen::ArrayXd::Zero(113);
+    Eigen::ArrayXd TORQUE = Eigen::ArrayXd::Zero(113);
+    for (size_t i = 0; i != SPEED.size(); ++i) {
+        SPEED(i) = (model.Params.CornerSpeed + i * (model.Params.MaxSpeed - model.Params.CornerSpeed) / (SPEED.size() - 1)) * 0.030 / M_PI;
+    }
+
+    // Series mode torque-speed curve
+    //get_torque_and_flux_linkage(model.msp, model.position, Ns, model.pp, torque, flux_cos, flux_sin, angle, current);
+    model.simulate(torque, flux_cos, flux_sin, angle, current);
+
+    double_t max_torque = torque.array().maxCoeff();
+
+    torque_speed_curve(torque,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE);
+
+    // Parallel mode torque-speed curve
+    current.array() *= 0.5;
+
+    //get_torque_and_flux_linkage(model.msp, model.position, Ns, model.pp, torque, flux_cos, flux_sin, angle, current);
+    model.simulate(torque, flux_cos, flux_sin, angle, current);
+
+    for (size_t i = 0; i != 3; ++i) {
+        flux_cos[i].array() *= 0.5;
+        flux_sin[i].array() *= 0.5;
+    }
+
+    torque_speed_curve(torque,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE);
+
+    // Output
+    Eigen::ArrayXd POWER = TORQUE * SPEED * M_PI / 30.0;
+    /*
+    for (size_t i = 0; i != SPEED.size(); ++i) {
+        std::cout << SPEED[i] << "," << (size_t)TORQUE[i] << "," << (size_t)POWER[i] << std::endl;
+    }
+    */
+
+    // Post Process
+    ObjectiveMap om;
+
+    om["power"] = POWER.minCoeff();
+    om["torque"] = max_torque;
+
+    double_t backiron_radius = model.dws.outer_radius();
+    double_t active_volume = model.Params.ls * pow(backiron_radius,2.0) * M_PI * 1e3;
+
+    om["active volume"] = active_volume;
+
+    double_t end_turn_length = 2 * model.dws.outer_radius() * sin(M_PI / 2.0 / model.Params.Np);
+    double_t total_length = end_turn_length + model.Params.ls;
+    double_t total_volume = total_length * M_PI * pow(backiron_radius,2.0) * 1e3;
+
+    om["total volume"] = total_volume;
+
+    //std::cout << min_power << "," << max_torque << "," << active_volume << "," << total_volume << std::endl;
+
+    return om;
 }
 
 std::function<ObjectiveMap(Particle)> global_objective() {

apps/Switched Winding Synchrel/main.cpp

 #include "designspace.h"
 #include "globalobjective.h"
+#include "model.h"
+
+#define SAVE_DIR "./"
 
 int main(int argc, char** argv) {
+    MPI_Init(NULL,NULL);
+    int rank{0};
+    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+
+    int size{0};
+    MPI_Comm_size(MPI_COMM_WORLD,&size);
+
     // Arguments
-    size_t swarm_size{41};
-    double_t objective_tolerance{1e-2};
-    size_t maximum_iterations{20};
+    size_t swarm_size{31};
+    double_t objective_tolerance{1e-2 * 22.0};
+    size_t maximum_iterations{3};
 
     CoordinateSpace ds = design_space();
-    std::cout << ds;
+
+    if (rank == 0) {
+        std::cout << ds;
+    }
 
     std::function<ObjectiveMap(Particle)> go = global_objective();
 
-    Swarm swarm{ds,go,swarm_size,objective_tolerance,maximum_iterations};
+    Swarm swarm{ds,go,swarm_size,objective_tolerance,maximum_iterations,std::greater<double_t>()};
+
+    for (auto &p : swarm.Particles) {
+        p.local_objective() = [](ObjectiveMap m) {
+            if (m["power"] < 55.0) {
+                return ((m["power"] - 55.0) / m["active volume"]);
+            } else {
+                return m["power"] / m["active volume"];
+            }
+        };
+    }
+
+    if(size > 1) {
+        swarm.mpi_run();
+    } else {
+        swarm.run();
+    }
 
-    std::cout << swarm;
+    if (rank == 0) {
+        std::cout << swarm;
+        std::cerr << "Reminder: Implement analytical bridge thickness model" << std::endl;
+    }
 
-    std::cerr << "Reminder: Don't forget to group limit at0-at5 and rt0-rt6" << std::endl;
+    MPI_Finalize();
 
     return 0;
 }
\ No newline at end of file

apps/Switched Winding Synchrel/model.cpp

+#include "model.h"
+
+void Model::create_sketch() {
+    srr.Poles = Params.Np;
+    srr.InnerRadius = Params.rri;
+    srr.OuterRadius = Params.rro;
+
+    srr.AngularBridgeThickness = Params.angular_bridge_thickness;
+    srr.RadialBridgeThickness = Params.radial_bridge_thickness;
+    srr.RadialThickness = Params.radial_thickness;
+    srr.AngularThickness = Params.angular_thickness;
+
+    srr.Convexify = true;
+    srr.AddBoundingBox = false;
+
+    dws.InnerRadius = Params.rsi;
+    dws.ToothFaceThickness = Params.tooth_face_thickness;
+    dws.SlotDepth = Params.slot_depth;
+    dws.BackironThickness = Params.backiron_thickness;
+
+    dws.SlotWidth = Params.slot_width;
+    dws.SlotOpening = Params.slot_opening;
+
+    dws.TotalTeeth = Params.Nt;
+    dws.ModeledTeeth = Params.Nt / Params.Np;
+
+    dws.Convexify = false;
+    dws.AddBoundingBox = true;
+
+    origin = sketch.new_element<Vertex>(0.0,0.0);
+    sketch.new_element<Fixation>(origin);
+
+    srr.add_to_sketch(sketch,origin);
+    dws.add_to_sketch(sketch,origin);
+
+    airgap = CylindricalAirgap<2>{sketch,origin,Params.rro,0.0,Params.rsi,0.0,360.0/Params.Np};
+
+    double_t tol = sketch.solve();
+    if (tol > 100e-3 * FLT_EPSILON) {
+        std::cerr << "Sketch was not solved within specified tolerance" << std::endl;
+    }
+    //EXPECT_LE(sketch.solve(), 100e-3 * FLT_EPSILON);
+    //EXPECT_TRUE(sketch.build());
+
+    bool success = sketch.build();
+    if (!success) {
+        std::cerr << "Sketch was not built successfully" << std::endl;
+    }
+
+    // Sketch selections for boundary conditions
+    airgap_motion_interface = airgap.interface(1);
+    periodic_boundaries = sketch.select_periodic_boundary_pairs(origin, 360.0 / Params.Np);
+    magnetic_insulation_boundary = dws.bounding_box_circular_arc();
+}
+
+void Model::create_mesh() {
+    mesh = Mesh{sketch};
+
+    mesh.boundary_constraint(airgap_motion_interface)->uniform_discretization(true);
+    mesh.add_mapped_boundary_pair(periodic_boundaries);
+
+    mesh.create();
+
+    if(!mesh.edges_are_valid()) {
+        std::cerr << "Initial mesh edges are not valid" << std::endl;
+    }
+
+    if(!mesh.edges_are_optimal()) {
+        std::cerr << "Initial mesh edges are not optimal" << std::endl;
+    }
+
+    mesh.MinimumElementQuality = 0.2;
+
+    if(!mesh.refine()) {
+        std::cerr << "Mesh was not successfully refined" << std::endl;
+    }
+
+
+    if(!mesh.edges_are_valid()) {
+        std::cerr << "Refined mesh edges are not valid" << std::endl;
+    }
+
+    if(!mesh.edges_are_valid()) {
+        std::cerr << "Refined mesh edges are not optimal" << std::endl;
+    }
+
+    // Mesh selections for material properties and forcing functions
+    double_t r = srr.inner_radius() * (1 + 1e-3);
+    double_t x = r * std::cos(M_PI/Params.Np);
+    double_t y = r * std::sin(M_PI/Params.Np);
+    rotor_iron = mesh.select_contour(Point{x,y});
+
+    r = dws.outer_radius() * (1 - 1e-3);
+    x = r * std::cos(M_PI/Params.Np);
+    y = r * std::sin(M_PI/Params.Np);
+    stator_iron = mesh.select_contour(Point{x,y});
+
+    slot_contours = dws.select_slot_contours(mesh);
+
+    r = (srr.OuterRadius + dws.InnerRadius) / 2.0;
+    airgap_contour = mesh.select_contour(Point{r*cos(M_PI / Params.Np),r*sin(M_PI / Params.Np)});
+}
+
+void Model::create_physics() {
+    // Setup FEA simulation
+    fem = std::make_shared<FiniteElementMesh<2, 2>>(mesh);
+    msp = std::make_shared<Magnetostatic>(fem);
+
+    // Set material properties
+    MaterialProperties electrical_steel = JFE_35JN210();
+    fem->region(rotor_iron)->material(electrical_steel);
+    fem->region(stator_iron)->material(electrical_steel);
+
+    // Current Density
+    msp->add_current_density([](FunctionArguments args) { return -args["J"] * cos(2 * M_PI * args["t"] - M_PI * 2 / 3 + args["a"]);}, {slot_contours[0],slot_contours[1]});
+    msp->add_current_density([](FunctionArguments args) { return +args["J"] * cos(2 * M_PI * args["t"] + M_PI * 0 / 3 + args["a"]);}, {slot_contours[2],slot_contours[3]});
+    msp->add_current_density([](FunctionArguments args) { return -args["J"] * cos(2 * M_PI * args["t"] + M_PI * 2 / 3 + args["a"]);}, {slot_contours[4],slot_contours[5]});
+
+    // Boundary Conditions
+    msp->add_magnetic_insulation({magnetic_insulation_boundary});
+
+    // Rotating Interface
+    position = msp->add_sliding_interface(airgap_motion_interface, true);
+    double_t dt = 2.0 / (4.0 * position->size());
+
+    msp->add_periodic_boundary(periodic_boundaries, true);
+}
+
+void Model::create_post_processor() {
+    pp.add<AirgapTorque>(msp,"Torque",fem->region(airgap_contour), Params.rro*2.0/3.0 + Params.rsi*1.0/3.0, Params.rro*1.0/3.0 + Params.rsi*2.0/3.0, 1.0, Params.Np * Params.ls);
+    pp.add<Integral>(msp,"FluxA",fem->region({slot_contours[0],slot_contours[1]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax * Params.ls, false);
+    pp.add<Integral>(msp,"FluxB",fem->region({slot_contours[2],slot_contours[3]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax* Params.ls, false);
+    pp.add<Integral>(msp,"FluxC",fem->region({slot_contours[4],slot_contours[5]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax* Params.ls, false);
+}
+
+void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current) {
+    torque.setZero();
+    for (size_t i = 0; i!=3; ++i) {
+        flux_cos[i].setZero();
+        flux_sin[i].setZero();
+    }
+
+    size_t Ns{16};
+    double_t t{0.0};
+    double_t dt = 2.0 / (4.0 * position->size());
+    size_t dp{position->size() / Ns};
+
+    std::vector<Eigen::VectorXd> previous(angle.size());
+    for (size_t is = 0; is != Ns; ++is) {
+        msp->assemble();
+
+#pragma omp parallel for
+        for (size_t i = 0; i < angle.size(); i++) {
+            auto solution = msp->new_solution();
+            if (is > 1)
+                solution->v = previous[i];
+
+            FunctionArguments fargs;
+            fargs.add_argument("t",t);
+            fargs.add_argument("a",angle(i));
+            fargs.add_argument("J",current(i));
+
+            msp->solve(solution, fargs);
+            previous[i] = solution->v;
+
+            //Verify equation is solved
+            //EXPECT_LE(solution->r.norm(), 1e-2 * (solution->J * solution->v).norm());
+
+            torque(i) += pp("Torque", solution->v) / Ns;
+            flux_cos[0](i) += pp("FluxA",solution->v) * (2.0 * cos(2.0 * M_PI * t) / Ns);
+            flux_cos[1](i) += pp("FluxB",solution->v) * (2.0 * cos(2.0 * M_PI * t) / Ns);
+            flux_cos[2](i) += pp("FluxC",solution->v) * (2.0 * cos(2.0 * M_PI * t) / Ns);
+
+            flux_sin[0](i) += pp("FluxA",solution->v) * (2.0 * sin(2.0 * M_PI * t) / Ns);
+            flux_sin[1](i) += pp("FluxB",solution->v) * (2.0 * sin(2.0 * M_PI * t) / Ns);
+            flux_sin[2](i) += pp("FluxC",solution->v) * (2.0 * sin(2.0 * M_PI * t) / Ns);
+        }
+
+        // Increment Position
+        t += dp * dt;
+        *position += dp;
+    }
+    *position += Ns * dp;
+}

apps/Switched Winding Synchrel/model.h

+#ifndef OE_APP_SWITCHED_WINDING_SYNCHREL_MODEL_H
+#define OE_APP_SWITCHED_WINDING_SYNCHREL_MODEL_H
+
+#include "Sketch.hpp"
+#include "Mesh.hpp"
+#include "Physics.hpp"
+#include "ModelTemplates.hpp"
+#include "Optimization.hpp"
+
+#include "parameters.h"
+
+class Model {
+public:
+    Model(Particle particle) : Params{particle} {
+        create_sketch();
+        create_mesh();
+        create_physics();
+        create_post_processor();
+    };
+
+    Parameters Params;
+
+    // Sketch
+    Sketch sketch;
+    SynchronousReluctanceRotor srr;
+    DistributedWindingStator dws;
+
+    std::shared_ptr<Vertex> origin;
+    CylindricalAirgap<2> airgap;
+
+    std::shared_ptr<CircularArc> airgap_motion_interface;
+    std::vector<PeriodicBoundaryPair> periodic_boundaries;
+    std::shared_ptr<Curve const> magnetic_insulation_boundary;
+
+    // Mesh
+    Mesh mesh;
+
+    std::shared_ptr<const Contour> rotor_iron;
+    std::shared_ptr<const Contour> stator_iron;
+    std::shared_ptr<const Contour> airgap_contour;
+    std::vector<std::shared_ptr<const Contour>> slot_contours;
+
+    // Physics
+    std::shared_ptr<FiniteElementMesh<2,2>> fem;
+    std::shared_ptr<Magnetostatic> msp;
+
+    std::shared_ptr<SlidingInterface> position;
+
+    // Post Processor
+    PostProcessorInterface pp;
+
+    void simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current);
+
+protected:
+    void create_sketch();
+
+    void create_mesh();
+
+    void create_physics();
+
+    void create_post_processor();
+};
+
+#endif //OE_APP_SWITCHED_WINDING_SYNCHREL_MODEL_H

apps/Switched Winding Synchrel/parameters.cpp

+#include "parameters.h"
+
+Parameters::Parameters(Particle p) {
+    std::cerr << "Reminder: Implement analytical bridge thickness model" << std::endl;
+
+    ls = p["ls"].Position;
+    rro = p["rro"].Position;
+
+    rsi = rro + 0.03*INCH; // airgap length 30mils
+
+    angular_bridge_thickness = {0.5e-3,0.5e-3,0.5e-3};
+    radial_bridge_thickness = {0.5e-3,0.5e-3,0.5e-3};
+
+    std::vector<double_t> rt(6,0.0);
+    for (size_t i = 0; i != 6; ++i) {
+        std::string key = "rt" + std::to_string(i);
+        rt[i] = 2 * rro * sin(M_PI / 2 / Np) * p[key].Position;
+    }
+    radial_thickness = rt;
+
+    std::vector<double_t> at(6,0.0);
+    for (size_t i = 0; i != 6; ++i) {
+        std::string key = "at" + std::to_string(i);
+        at[i] = 180.0 / Np * p[key].Position;
+    }
+    angular_thickness = at;
+
+    slot_width = p["sw"].Position * (2 * M_PI * rro) / (Nt);
+
+    slot_depth = (p["sd"].Position / p["sw"].Position) * (Imax * Nt) / (2 * M_PI * rro * Jmax * SlotFill);
+
+    backiron_thickness = p["bi"].Position * (2 * M_PI * (rro + slot_depth - rro * (1 - p["sw"].Position))) / (2 * Np);
+}
\ No newline at end of file

apps/Switched Winding Synchrel/parameters.h

+#ifndef OE_APP_SWITCHED_WINDING_SYNCHREL_PARAMETERS_H
+#define OE_APP_SWITCHED_WINDING_SYNCHREL_PARAMETERS_H
+
+#include "Optimization.hpp"
+
+class Parameters {
+public:
+    Parameters(Particle p);
+
+    double_t Vdc{600.0}; // dc link voltage
+    double_t Imax{400.0*325.0/600.0}; // maximum phase current
+    double_t Jmax{30.0e6}; // maximum current density
+    double_t SlotFill{0.5}; // slot fill factor
+
+    double_t MaxSpeed{14.0e3*M_PI/30.0}; // rad/s
+    double_t CornerSpeed{2.8e3*M_PI/30.0}; // rad/s
+
+    size_t Np{8};
+    size_t Nt{48};
+
+    double_t ls;
+    double_t rri{25.4e-3};
+    double_t rro;
+    double_t rsi;
+
+    std::vector<double_t> angular_bridge_thickness;
+    std::vector<double_t> radial_bridge_thickness;
+    std::vector<double_t> radial_thickness;
+    std::vector<double_t> angular_thickness;
+
+    double_t tooth_face_thickness{2*0.35e-3};
+    double_t slot_depth;
+    double_t backiron_thickness;
+
+    double_t slot_width;
+    double_t slot_opening{1.74e-3};
+protected:
+    double_t INCH{25.4e-3};
+};
+
+#endif //OERSTED_PARAMETERS_H

apps/Switched Winding Synchrel/torquespeedcurve.cpp

+#include "torquespeedcurve.h"
+
+void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, double_t DC_LINK, double_t Np, Eigen::ArrayXd &SPEED, Eigen::ArrayXd &TORQUE) {
+    size_t Nr = (size_t)torque.rows();
+    size_t Nc = (size_t)torque.cols();
+
+    Eigen::MatrixXd phase_flux = Eigen::MatrixXd::Zero(Nr,Nc);
+    Eigen::MatrixXd max_speed = Eigen::MatrixXd::Zero(Nr,Nc);
+
+    for (size_t i = 0; i != torque.rows(); ++i) {
+        for (size_t j = 0; j != torque.cols(); ++j) {
+            for (size_t k = 0; k != 3; ++k) {
+                phase_flux(i,j) += std::hypot(flux_cos[k](i,j), flux_sin[k](i,j));
+            }