Commit 5d7843f2 authored by JasonPries's avatar JasonPries
Browse files

Nightly Commit

Initial implementation and testing of periodic/antiperiodic boundary conditions for magnetostatic solvers
parent d093a288
......@@ -25,9 +25,59 @@ protected:
std::vector<size_t> Nodes;
};
/*
template <size_t Dimension>
class BoundaryPair {
class BoundaryPair {};
template <>
class BoundaryPair<2> {
public:
BoundaryPair(std::shared_ptr<Boundary<2>> first, std::shared_ptr<Boundary<2>> second) : First{first}, Second{second} {};
std::shared_ptr<Boundary<2>> const &first() { return First; };
std::shared_ptr<Boundary<2>> const &second() { return Second; };
protected:
std::shared_ptr<Boundary<2>> First;
std::shared_ptr<Boundary<2>> Second;
};
*/
class VariableMap {
public:
VariableMap(size_t first, size_t second, double value) : First{first}, Second{second}, Value{value} {};
size_t first() const { return First; };
size_t second() const { return Second; };
double value() const { return Value; };
protected:
size_t First;
size_t Second;
double Value;
};
class BoundaryMap {
public:
BoundaryMap(std::vector<VariableMap> map) : Map{map} {};
size_t size() const { return Map.size(); };
VariableMap const &map(size_t i) const { return Map[i]; };
std::vector<VariableMap> const &map() const { return Map; };
size_t first(size_t i) const { return Map[i].first(); };
size_t second(size_t i) const { return Map[i].second(); };
double value(size_t i) const { return Map[i].value(); };
protected:
std::vector<VariableMap> Map;
};
#endif //OERSTED_BOUNDARY_H
......@@ -16,25 +16,27 @@
class BoundaryCondition {
public:
virtual ~BoundaryCondition(){};
virtual ~BoundaryCondition() {};
virtual void apply(std::vector<Eigen::Triplet<double>> &triplets) const = 0;
virtual void reduce(std::set<size_t,std::greater<size_t>> &index) const = 0;
virtual void reduce(std::set<size_t, std::less<size_t>> &index) const = 0;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet : public BoundaryCondition{};
class ZeroDirichlet : public BoundaryCondition {
};
template<size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
class ZeroDirichlet<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
public:
ZeroDirichlet(std::vector<std::shared_ptr<Boundary<2>>> boundaries) : Boundaries{boundaries} {};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};
void reduce(std::set<size_t,std::greater<size_t>> &index) const override {
for (auto const& b : Boundaries) {
for (auto const& i : b->nodes()) {
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (auto const &b : Boundaries) {
for (auto const &i : b->nodes()) {
index.insert(i);
}
}
......@@ -46,37 +48,34 @@ protected:
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition : public BoundaryCondition{};
class PeriodicBoundaryCondition : public BoundaryCondition {
};
template<size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
public:
PeriodicBoundaryCondition(std::vector<std::shared_ptr<BoundaryPair<2>>> boundary_pairs, bool antiperiodic) : BoundaryPairs{boundary_pairs}, Antiperiodic{antiperiodic} {};
PeriodicBoundaryCondition(std::vector<BoundaryMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
double sign{1.0 - 2.0 * Antiperiodic};
// TODO:
/*
for(size_t i = 0;i!=FirstNodes.size();++i){
triplets[FirstNodes[i]].row() = SecondNodes[i];
triplest[FirstNodes[i]].value() = sign;
}
*/
for (BoundaryMap const &b : Map) {
for (VariableMap const &v : b.map()) {
triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), sign * v.value()));
}
}
}
void reduce(std::set<size_t, std::greater<size_t>> &index) const override {
/*
for (auto const &bp : BoundaryPairs) {
for(auto const &ij : bp->nodes()) {
index.insert(ij.second);
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (BoundaryMap const &b : Map) {
for (VariableMap const &v : b.map()) {
index.insert(v.second());
}
}
*/
}
protected:
std::vector<std::shared_ptr<BoundaryPair<2>>> BoundaryPairs;
std::vector<BoundaryMap> Map;
bool Antiperiodic;
};
......
......@@ -61,9 +61,9 @@ public:
*/
PhysicsData(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
ElementWeights{d.triangles().size()},
Derivatives{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes},
Basis{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes} {};
ElementWeights{d.triangles().size()},
Derivatives{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes},
Basis{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes} {};
FiniteElementMesh<Dimension, ElementOrder> const &domain() const { return Domain; };
......@@ -120,8 +120,8 @@ public:
Fy = -Derivatives.dx(i).transpose() * v;
// Calculate polarization
for(auto &dr : Domain.regions()) {
dr->material().h_from_b(dr->elements(),Fx,Fy,dFxdx,dFydy,dFxdy);
for (auto &dr : Domain.regions()) {
dr->material().h_from_b(dr->elements(), Fx, Fy, dFxdx, dFydy, dFxdy);
}
// Calculate residual
......@@ -149,7 +149,7 @@ public:
void calculate_forcing(Eigen::VectorXd &f) {
f.setZero();
for(size_t i = 0; i!=ForcingCondtions.size();++i) {
for (size_t i = 0; i != ForcingCondtions.size(); ++i) {
f += ForcingCondtions[i]->operator()(Time);
}
}
......@@ -166,52 +166,60 @@ public:
BoundaryConditions.push_back(std::make_unique<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
}
void add_periodic_boundary(std::vector<BoundaryMap> map, bool antiperiodic) {
BoundaryConditions.push_back(std::make_unique<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic));
}
void apply_conditions() {
// TODO: This is a general method that should be in a superclass
std::vector<Eigen::Triplet<double>> triplets;
triplets.reserve(Domain.nodes().size());
for(size_t i = 0; i!=Domain.nodes().size();++i) {
triplets.push_back(Eigen::Triplet<double>(i,i,1.0));
for (size_t i = 0; i != Domain.nodes().size(); ++i) {
triplets.push_back(Eigen::Triplet<double>(i, i, 1.0));
}
for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
BoundaryConditions[i]->apply(triplets);
}
std::set<size_t, std::greater<size_t>> index;
std::set<size_t, std::less<size_t>> index;
for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
BoundaryConditions[i]->reduce(index);
}
// TODO: Rewrite following section:
// (Need to handle arbitrary number of entries in the b.c. "permutation" matrix)
// First, search for triplets with row() == reduce and erase()
// Then, search for skipped rows and cumulatively renumber
// BEGIN
size_t rows{Domain.nodes().size()};
size_t cols{Domain.nodes().size()};
for (auto &i : index) { // index is sorted from greatest to least, so triplets are erased in reverse order
auto iter = (triplets.begin() + i);
// TODO: Does this handle arbitrary number of entries in the b.c. "permutation" matrix correctly?
// TODO: BEGIN
// Sort triplet list by rows
auto row_lt = [](Eigen::Triplet<double> const &x, Eigen::Triplet<double> const &y) { return x.row() < y.row(); };
std::sort(triplets.begin(), triplets.end(), row_lt);
size_t num_erased{0};
auto iter{triplets.begin()};
for (size_t const &i : index) {
// Renumber triplets based on number of rows previously erased
while (iter != triplets.end() && iter->row() < i) {
*iter = Eigen::Triplet<double>(iter->row() - num_erased, iter->col(), iter->value());
++iter;
}
if(iter->row() == i) {
triplets.erase(iter);
// Erase all rows associated with the
while (iter != triplets.end() && iter->row() == i) {
iter = triplets.erase(iter);
}
--rows;
++num_erased;
}
for (size_t i = 0; i != triplets.size(); ++i) {
triplets[i] = Eigen::Triplet<double>(i, triplets[i].col(), triplets[i].value());
}
// END
size_t cols{Domain.nodes().size()};
size_t rows{cols - num_erased};
// TODO: END
Eigen::SparseMatrix<double> bc_matrix(rows,cols);
bc_matrix.setFromTriplets(triplets.begin(),triplets.end());
Eigen::SparseMatrix<double> bc_matrix(rows, cols);
bc_matrix.setFromTriplets(triplets.begin(), triplets.end());
Basis.transform(bc_matrix);
Derivatives.transform(bc_matrix);
Unknowns = rows;
};
......
......@@ -451,7 +451,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
*/
}
class TwoRegionHex : public ::testing::Test {
class TwoRegionHexagon : public ::testing::Test {
public:
virtual void SetUp() {
node.push_back(XY{0.0, 0.0});
......@@ -568,7 +568,7 @@ public:
}
};
TEST_F(TwoRegionHex, magnetostatic_forcing_air) {
TEST_F(TwoRegionHexagon, magnetostatic_forcing_air) {
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
ms.time(0.0);
......@@ -636,7 +636,7 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_air) {
}
// Fundamental solution is 0.25*r^2 + C in forced region (0)
EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, FLT_EPSILON);
EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * FLT_EPSILON);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
}
......@@ -656,7 +656,7 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_air) {
test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
}
TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
TEST_F(TwoRegionHexagon, magnetostatic_forcing_1000) {
femesh.region(1)->material(Linear1000);
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
......@@ -698,7 +698,7 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
}
// Fundamental solution is 0.25*r^2 + C in forced region (0)
EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, FLT_EPSILON);
EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * FLT_EPSILON);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
}
......@@ -724,4 +724,219 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
}
test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
}
class TwoRegionQuarterOctagon : public ::testing::Test {
public:
virtual void SetUp() {
// Inner Region
node.push_back(XY{0.0, 0.0});
double r{1.0};
for (size_t i = 0; i != 3; ++i) {
double a = (i * M_PI) / 4.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
inner_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3});
bnd.push_back(inner_boundary);
tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
tri.push_back(Triangle<1>(1, 0_zu, 2_zu, 3_zu));
inner_region = std::make_shared<Region<2>>(std::vector<size_t>{0, 1});
reg.push_back(inner_region);
// Outer Region
r = 2.0;
for (size_t i = 0; i != 3; ++i) {
double a = (i * M_PI) / 4.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
outer_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{4, 5, 6});
bnd.push_back(outer_boundary);
tri.push_back(Triangle<1>(2, 1_zu, 4_zu, 2_zu));
tri.push_back(Triangle<1>(3, 4_zu, 5_zu, 2_zu));
tri.push_back(Triangle<1>(4, 2_zu, 5_zu, 6_zu));
tri.push_back(Triangle<1>(5, 2_zu, 6_zu, 3_zu));
outer_region = std::make_shared<Region<2>>(std::vector<size_t>{2, 3, 4, 5});
reg.push_back(outer_region);
// Periodic Boundary
variable_map.push_back(VariableMap{1, 3, 1.0});
variable_map.push_back(VariableMap{4, 6, 1.0});
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<std::shared_ptr<Region<2>>> reg;
std::vector<std::shared_ptr<Boundary<2>>> bnd;
std::shared_ptr<Region<2>> inner_region;
std::shared_ptr<Region<2>> outer_region;
std::shared_ptr<Boundary<2>> inner_boundary;
std::shared_ptr<Boundary<2>> outer_boundary;
std::vector<VariableMap> variable_map;
FiniteElementMesh<2, 1> femesh;
};
TEST_F(TwoRegionQuarterOctagon, magnetostatic_periodic) {
// Initialize physics
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
ms.time(0.0);
// Add conditions
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, {inner_region});
ms.add_magnetic_insulation({outer_boundary});
ms.add_periodic_boundary({variable_map}, false);
ms.build_matrices();
ms.apply_conditions();
// Initialize matrices
auto J = ms.init_unknown_matrix();
auto v = ms.init_unknown_vector();
auto r = ms.init_unknown_vector();
auto f = ms.init_unknown_vector();
auto Fx = ms.init_element_array();
auto Fy = ms.init_element_array();
auto dFxdx = ms.init_element_array();
auto dFydy = ms.init_element_array();
auto dFxdy = ms.init_element_array();
v.setZero();
ms.calculate_forcing(f);
double nodal_current{M_SQRT2 / 6.0};
for (size_t i = 0; i != 3; ++i) {
EXPECT_NEAR(f(i), nodal_current, FLT_EPSILON);
}
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
ldlt.compute(J);
ASSERT_EQ(ldlt.info(), Eigen::Success);
v -= ldlt.solve(r);
// Verify equation is solved
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
for (size_t i = 0; i != r.size(); ++i) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
}
// Fundamental solution is 0.25*r^2 + C in forced region (0)
EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * 0.1);
for (size_t i = 1; i != 3; ++i) {
EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
}
// Test flux-density values
Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
Eigen::ArrayXd Bmag(By.size());
Eigen::ArrayXd Bang(By.size());
for (size_t i = 0; i != By.size(); ++i) {
Bmag(i) = hypot(Bx(i), By(i));
Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
}
std::vector<double> b_mag_expected(6, 3.86994e-7);
std::vector<double> b_ang_expected{112.5, 157.5, 112.5, 112.5, 157.5, 157.5};
for(size_t i = 0;i!=6;++i) {
EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-12);
EXPECT_NEAR(Bang[i], b_ang_expected[i], 1e-1);
}
}
TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) {
// Initialize physics
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
ms.time(0.0);
// Add conditions
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, {inner_region});
ms.add_magnetic_insulation({outer_boundary});
variable_map.push_back(VariableMap(0, 0, -1.0));
ms.add_periodic_boundary({variable_map}, true);
ms.build_matrices();
ms.apply_conditions();
// Initialize matrices
auto J = ms.init_unknown_matrix();
auto v = ms.init_unknown_vector();
auto r = ms.init_unknown_vector();
auto f = ms.init_unknown_vector();
auto Fx = ms.init_element_array();
auto Fy = ms.init_element_array();
auto dFxdx = ms.init_element_array();
auto dFydy = ms.init_element_array();
auto dFxdy = ms.init_element_array();
v.setZero();
ms.calculate_forcing(f);
EXPECT_NEAR(f(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(f(1), M_SQRT2 / 6.0, FLT_EPSILON);
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
ldlt.compute(J);
ASSERT_EQ(ldlt.info(), Eigen::Success);
v -= ldlt.solve(r);
// Verify equation is solved
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
for (size_t i = 0; i != r.size(); ++i) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
}
std::vector<double> v_expected{0.0,6.60365e-8};
for (size_t i = 0; i != 2; ++i) {
EXPECT_NEAR(v[i], v_expected[i], 1e-13);
}
// Test flux-density values
Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
Eigen::ArrayXd Bmag(By.size());
Eigen::ArrayXd Bang(By.size());
for (size_t i = 0; i != By.size(); ++i) {
Bmag(i) = hypot(Bx(i), By(i));
Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
}
std::vector<double> b_mag_expected{9.33897e-8, 9.33897e-8, 9.33897e-8, 7.14774e-8, 7.14774e-8, 9.33897e-8};
std::vector<double> b_ang_expected{0.0, -90, 0.0, 112.5, 157.5, -90.0};
for(size_t i = 0; i!= 6;++i) {
EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-13);
EXPECT_NEAR(Bang[i], b_ang_expected[i], FLT_EPSILON);
}
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment