Merge branch 'devel' into 'alpha'

Merge development branch as alpha release

See merge request p7k/Oersted!7

.idea/codeStyleSettings.xml

@@ -30,6 +30,6 @@
         </Objective-C-extensions>
       </value>
     </option>
-    <option name="PREFERRED_PROJECT_CODE_STYLE" value="Default (1)" />
+    <option name="PREFERRED_PROJECT_CODE_STYLE" value="Default" />
   </component>
 </project>
\ No newline at end of file

.idea/dictionaries/jpries.xml

 <component name="ProjectDictionaryState">
   <dictionary name="jpries">
     <words>
+      <w>discretization</w>
       <w>eigen</w>
+      <w>linearize</w>
       <w>supremum</w>
       <w>tangency</w>
       <w>verticies</w>

CMakeLists.txt

@@ -2,23 +2,37 @@ cmake_minimum_required(VERSION 3.2)
 
 project(Oersted)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++14")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage")
-#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
-
-find_package(Boost REQUIRED COMPONENTS system filesystem)
-include_directories(${Boost_INCLUDE_DIR})
+find_package(OpenMP REQUIRED)
+find_package(Boost REQUIRED)
+find_package(MPI REQUIRED)
 
+include_directories(SYSTEM ${MPI_CXX_INCLUDE_PATH})
 include_directories(./lib/)
 include_directories(./lib/Eigen/)
-include_directories(./lib/GoogleTest/googletest/include/gtest)
+include_directories(./lib/GoogleTest/googletest/include/gtest/)
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++14")
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
 
 add_subdirectory(./lib/)
 
 include_directories(./src/)
+include_directories(./src/IO/include/)
+include_directories(./src/Geometry/include/)
 include_directories(./src/Sketch/include/)
 include_directories(./src/Mesh/include/)
+include_directories(./src/Quadrature/include/)
+include_directories(./src/Elements/include/)
+include_directories(./src/Matrix/include/)
+include_directories(./src/Materials/include/)
+include_directories(./src/Optimization/include/)
+include_directories(./src/Physics/include/)
+include_directories(./src/Model\ Templates/include/)
 
 add_subdirectory(./src/)
 
-add_subdirectory(./test/)
\ No newline at end of file
+add_subdirectory(./test/)
+
+add_subdirectory(./apps/)
\ No newline at end of file

apps/CMakeLists.txt

+project(Oersted_Apps)
+
+set(OERSTED_LIBRARIES
+        io
+        geometry
+        sketch
+        elements
+        mesh
+        materials
+        optimization
+        quadrature
+        physics
+        model_templates
+        stdc++fs
+        ${MPI_CXX_LIBRARIES}
+        )
+
+add_subdirectory(./Switched\ Winding\ Synchrel)
\ No newline at end of file

apps/Switched Winding Synchrel/CMakeLists.txt

+PROJECT(Switched\ Winding\ Synchrel)
+
+set(SOURCE_FILES
+        main.cpp
+        designspace.h designspace.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 != 6; ++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 != 6; ++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 != 6; ++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;
+
+    double_t inch{25.4e-3};
+
+    d["Jmax"] = Coordinate{30e6,30e6,30e6};
+
+    d["ls"] = Coordinate{1*inch,3.5*inch,6*inch};
+    d["rro"] = Coordinate{2*inch,3.25*inch,4.5*inch};
+
+    for (size_t i = 0; i != 6; ++i) {
+        std::string key = "at" + std::to_string(i);
+        d[key] = Coordinate{0.07,0.14,0.21};
+    }
+
+    for (size_t i = 0; i != 6; ++i) {
+        std::string key = "rt" + std::to_string(i);
+        d[key] = Coordinate{0.07,0.14,0.21};
+    }
+
+    d["sw"] = Coordinate{0.2,0.4,0.6};
+    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;
+}

apps/Switched Winding Synchrel/designspace.h

+#ifndef OE_APP_SWITCHED_WINDING_SYNCHREL_DESIGNSPACE_H
+#define OE_APP_SWITCHED_WINDING_SYNCHREL_DESIGNSPACE_H
+
+#include "Optimization.hpp"
+
+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 particle) {
+    Model model{particle};
+
+    // Setup OpenMP
+    #ifdef _OPENMP
+        omp_set_num_threads(omp_get_max_threads());
+    #endif
+
+    size_t Na{11}, Nj{10};
+    Eigen::MatrixXd angle = Eigen::MatrixXd::Zero(Na,Nj);
+    Eigen::MatrixXd current = Eigen::MatrixXd::Zero(Na,Nj);
+    Eigen::MatrixXd torque = Eigen::MatrixXd::Zero(Na,Nj);
+    Eigen::MatrixXd ripple = 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::array<Eigen::MatrixXd,4> losses;
+    for (size_t i = 0; i != losses.size(); ++i) {
+        losses[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);
+    Eigen::ArrayXd RIPPLE = Eigen::ArrayXd::Zero(113);
+    Eigen::ArrayXd LOSSES = 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, ripple, flux_cos, flux_sin, angle, current, losses);
+    double_t max_torque = torque.array().maxCoeff();
+
+    torque_speed_curve(torque,ripple,losses,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE,RIPPLE,LOSSES);
+
+    // 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, ripple, flux_cos, flux_sin, angle, current, losses);
+
+    for (size_t i = 0; i != 3; ++i) {
+        flux_cos[i].array() *= 0.5;
+        flux_sin[i].array() *= 0.5;
+    }
+
+    torque_speed_curve(torque,ripple,losses,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE,RIPPLE,LOSSES);
+
+    // Output
+    Eigen::ArrayXd POWER = TORQUE * SPEED * M_PI / 30.0;
+
+    // Post Process
+    ObjectiveMap om;
+
+    om["power"] = POWER.minCoeff();
+    om["torque"] = max_torque;
+    om["ripple"] = RIPPLE.maxCoeff();
+
+    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;
+
+    om["losses"] = LOSSES.maxCoeff();
+
+    std::cout << om["torque"] << "," << om["ripple"] << "," << om["power"] << "," << om["losses"] << std::endl;
+
+    return om;
+}
+
+std::vector<std::string> global_objective_keys() {
+    return {"power","torque","ripple","active volume","total volume","losses"};
+}
+
+std::function<ObjectiveMap(Particle)> global_objective() {
+    return evaluate;
+}
\ No newline at end of file

apps/Switched Winding Synchrel/globalobjective.h

+#ifndef OE_APP_SWITCHED_WINDING_SYNCHREL_GLOBALOBJECTIVE_H
+#define OE_APP_SWITCHED_WINDING_SYNCHREL_GLOBALOBJECTIVE_H
+
+#include "Optimization.hpp"
+
+std::function<ObjectiveMap(Particle)> global_objective();
+
+std::vector<std::string> global_objective_keys();
+
+ObjectiveMap evaluate(Particle p);
+
+#endif // OE_APP_SWITCHED_WINDING_SYNCHREL_GLOBALOBJECTIVE_H
\ No newline at end of file

apps/Switched Winding Synchrel/main.cpp

+#include "designspace.h"
+#include "globalobjective.h"
+#include "model.h"
+#include "IO.hpp"
+
+#define SAVE_DIR save_directory() + "/Oersted/output/"
+
+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 N{11};
+    size_t swarm_size{((N - 2) * (N - 1)) / 2};
+    double_t objective_tolerance{1e-2};
+    size_t maximum_iterations{200};
+
+    CoordinateSpace ds = design_space();
+
+    if (rank == 0) {
+        std::cout << ds;
+    }
+
+    std::function<ObjectiveMap(Particle)> go = global_objective();
+
+    Swarm swarm{ds,go,swarm_size,objective_tolerance,maximum_iterations,global_objective_keys(),std::greater<double_t>()};
+
+    if (rank == 0) {
+        swarm.set_output_file(SAVE_DIR,"sws");
+    }
+
+    double_t dw{1.0/(N-1.0)};
+    double_t wp{0.0},wl{0.0},wr{0.0};
+    for (size_t i = 0; i < swarm_size; ++i) {
+        wp += dw;
+        if (wp + wl > 1.0 - dw + FLT_EPSILON) {
+            wp = dw;
+            wl += dw;
+        }
+        wr = 1.0 - wp - wl;
+
+        if (rank == 0) {
+        std::cout << wp << "," << wl << "," << wr << std::endl;
+        }
+
+        swarm.Particles[i].local_objective() = [wp,wl,wr](ObjectiveMap m) {
+            double_t wpd{wp / 9.8};
+            double_t wls{wl * 5.5};
+            double_t wrl{wr * 0.1};
+            if (m["power"] < 55.0) {
+                double_t value{0.0};
+                value += (std::min(0.0,m["power"] - 55.0)) / 55.0;
+                return value;
+            } else {
+                double_t value{0.0};
+                value += wpd * m["power"] / m["total volume"];
+                value += wls / m["losses"];
+                value += wrl * (m["torque"] / (m["torque"] - m["ripple"]));
+                return value;
+            }
+        };
+    }
+
+    if(size > 1) {
+        swarm.mpi_run();
+    } else {
+        swarm.run();
+    }
+
+    if (rank == 0) {
+        std::cout << swarm;
+    }
+
+    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);
+
+    double_t winding_length = Params.ls + M_PI * dws.outer_radius() * sin(M_PI / 2 / Params.Np); // including end turn length estimate
+    double_t resistivity = 1.68e-8 * (1 + 0.00404*(150 - 20));
+    pp.add<Integral>(msp,"ConductionLosses",fem->region({slot_contours[0]}), resistivity * Params.Nt * winding_length * 1e-3, false);
+}
+
+void Model::simulate(Eigen::MatrixXd &torque, Eigen::MatrixXd &ripple, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current, std::array<Eigen::MatrixXd,4> &losses) {
+    torque.setZero();
+
+    for (size_t i = 0; i != 3; ++i) {
+        flux_cos[i].setZero();
+        flux_sin[i].setZero();
+    }
+
+    for (size_t i = 0; i != 3; ++i) {
+        losses[i].setZero();
+    }
+
+    Eigen::MatrixXd torque_min = +FLT_MAX*Eigen::MatrixXd::Ones(torque.rows(),torque.cols());
+    Eigen::MatrixXd torque_max = -FLT_MAX*Eigen::MatrixXd::Ones(torque.rows(),torque.cols());
+
+    size_t Ns{32};
+    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());
+    std::vector<std::array<Eigen::VectorXd,3>> Bx_cos(angle.size());
+    std::vector<std::array<Eigen::VectorXd,3>> By_cos(angle.size());
+
+    std::vector<std::array<Eigen::VectorXd,3>> Bx_sin(angle.size());
+    std::vector<std::array<Eigen::VectorXd,3>> By_sin(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