Merge pull request #36 from OpenSEMBA/unshielded-wires

Multipolar expansion comparison with BEM

Author: lmdiazangulo

test/DriverTest.cpp

@@ -597,7 +597,7 @@ TEST_F(DriverTest, lansink2024_two_wires_using_multipolar_expansion)
 	EXPECT_NEAR(0.0, relError(expectedL11, computedL11), rTol);
 }
 
-TEST_F(DriverTest, lansink2024_two_wires_shifted_and_centered)
+TEST_F(DriverTest, lansink2024_fdtd_cell_shifted_and_centered)
 {
 	// From:
 	// Rotgerink, J.L. et al. (2024, September).
@@ -721,6 +721,10 @@ TEST_F(DriverTest, lansink2024_single_wire_multipolar_in_cell_parameters)
 		EXPECT_NEAR(0.0, relError(expectedL00, computedL00), rTol);
 	}
 
+	// Check that multipolar expansion for the bare wire produces a 1 V at the boundary.
+	auto a0 = inCell.magnetic.at(0).ab[0].first;
+	auto Va = a0 / (2 * M_PI) * log(1.0 / 1e-3);
+	EXPECT_NEAR(1.0, Va, 1e-3);
 }
 
 TEST_F(DriverTest, getCFromGeneralizedC_two_wires_open)

test/ElectrostaticSolverTest.cpp

@@ -791,9 +791,97 @@ TEST_F(ElectrostaticSolverTest, lansink2024_single_wire_L00_with_floating)
 
 	exportSolution(s, getTestCaseName());
 
-	auto expectedL00 = 320e-9; // Table 3 result has a mistake. This is the correct value.
+	// Table 3 result has a mistake. 
+	// This is the correct value obtained through personal communication.
+	auto expectedL00 = 320e-9;
 
 	// 
 	double rTol = 0.04;
 	EXPECT_NEAR(0.0, relError(expectedL00, computedL00), rTol);
+}
+
+TEST_F(ElectrostaticSolverTest, lansink2024_small_one_centered_bem_comparison)
+{
+	// Case from:
+	// Rotgerink, J.L. et al. (2024, September).
+	// Numerical Computation of In - cell Parameters for Multiwire Formalism in FDTD.
+	// In 2024 International Symposium on Electromagnetic Compatibility
+	// EMC Europe(pp. 334 - 339). IEEE.
+
+	// SMALL WIRE IS CENTERED AT ORIGIN.
+	// This test compares multipolar expansion results between Tulip and BEM.
+
+	const std::string CASE{ "lansink2024_small_one_centered" };
+	const std::string fn{ casesFolder() + CASE + "/" + CASE + ".pulmtln.in.json" };
+	auto model{ Parser{fn}.readModel() };
+
+	auto fp = Driver::loadFromFile(fn).getFloatingPotentials().electric;
+
+	SolverInputs p;
+	p.dirichletBoundaries = {
+		{
+			{1, fp(0,0)}, // Conductor 0,
+			{2, fp(0,1)}, // Conductor 1,
+		}
+	};
+	p.openBoundaries = { 3 };
+
+	ElectrostaticSolver s{ *model.getMesh(), p };
+	s.Solve();
+	
+
+	int order = 2;
+
+	// Tulip coefficients.
+	multipolarCoefficients ab(order + 1);
+	mfem::Vector origin({ 0.0, 0.0 });
+	for (int n = 0; n < order + 1; n++) {
+		ab[n] = {
+			s.getChargeMomentComponent(n, 0, origin),
+			s.getChargeMomentComponent(n, 1, origin)
+		};
+	}
+
+	// BEM coefficients.
+	multipolarCoefficients abBEM(order + 1);
+	abBEM[0] = { 1.0, 0.0 };
+	abBEM[1] = { 0.005, 0.0 };
+	abBEM[2] = { 8.7819e-5, 0 };
+
+	// Checks normalizing by total charge.
+	const double aTol = 1e-4;
+	auto a0 = ab[0].first;
+	for (int n = 0; n < order + 1; n++) {
+		EXPECT_NEAR(ab[n].first / a0,  abBEM[n].first,  aTol) << "for n=" << n;
+		EXPECT_NEAR(ab[n].second / a0, abBEM[n].second, aTol) << "for n=" << n;
+	}
+
+	// For debugging.
+	// Export FEM solution.
+	exportSolution(s, getTestCaseName());
+
+	// Export multipolar expansion projection.
+	auto multipolarSolution = s.getSolution();
+	std::function<double(const Vector&)> f =
+		std::bind(&multipolarExpansion, std::placeholders::_1, ab, origin);
+	FunctionCoefficient fc(f);
+	multipolarSolution.phi->ProjectCoefficient(fc);
+	s.setSolution(multipolarSolution);
+	exportSolution(s, getTestCaseName()+"_from_multipolar_expansion_in_origin");
+
+	// Computes and exports multipolar expansion with center of charge.
+	{
+		ElectrostaticSolver s{ *model.getMesh(), p };
+		s.Solve();
+		auto centerOfCharge = s.getCenterOfCharge();
+		auto mCoeff = s.getMultipolarCoefficients(order);
+		std::function<double(const Vector&)> f =
+			std::bind(&multipolarExpansion, std::placeholders::_1, mCoeff, centerOfCharge);
+		FunctionCoefficient fc(f);
+		auto multipolarSolution = s.getSolution();
+		multipolarSolution.phi->ProjectCoefficient(fc);
+		s.setSolution(multipolarSolution);
+		exportSolution(s, getTestCaseName() + "_from_multipolar_expansion_in_center_of_charge");
+	}
+
 }

testData/lansink2024_small_one_centered/lansink2024_small_one_centered.FCStd

@@ -0,0 +1,24 @@
+{
+   "geometries": [
+      {
+         "geometry": "Vacuum_0",
+         "area": 1066.539159843431
+      },
+      {
+         "geometry": "Vacuum_1",
+         "area": 67027.4826321312
+      },
+      {
+         "geometry": "Conductor_0",
+         "area": 3.141592653589793
+      },
+      {
+         "geometry": "Conductor_1",
+         "area": 314.1592653589793
+      },
+      {
+         "geometry": "OpenBoundary_0",
+         "area": 68411.3226499872
+      }
+   ]
+}