Merge pull request #41 from OpenSEMBA/numerical-comparisons

Numerical comparisons

Author: lmdiazangulo

.gitignore

@@ -68,3 +68,4 @@ __pycache__/
 /testData/lansink2024_fdtd_cell/calculations.py
 /testData/lansink2024/inCellPotentials.out.json
 /XMLGoogleTestResults.xml
+*.inCellPotentials.out.json

external/step2gmsh

@@ -19,6 +19,13 @@ struct Box {
 			(point[1] >= min[1] && point[1] <= max[1]);
 	}
 
+	void displace(const std::array<double,2>& point) {
+		for (int i = 0; i < 2; ++i) {
+			min[i] += point[i];
+			max[i] += point[i];
+		}
+	}
+
 	bool operator==(const Box& rhs) const
 	{
 		return min == rhs.min && max == rhs.max;

src/Model.h

@@ -271,7 +271,7 @@ double InCellPotentials::getInductanceOnBox(int i, int j, const Box& cellBox) co
     double Ij = magnetic.at(j).ab[0].first;
 
     double avAj = getAveragePotential(magnetic.at(j), innerRegionBox, cellBox);
-    double AiWhenPrescribedAj = electric.at(j).conductorPotentials.at(i);
+    double AiWhenPrescribedAj = magnetic.at(j).conductorPotentials.at(i);
     avAj = -avAj + AiWhenPrescribedAj;
 
     return avAj / Ij * MU0_SI;

src/Results.cpp

@@ -12,9 +12,9 @@ int main(int argc, char* argv[])
 	std::string inputFilename;
 
 	po::options_description desc(
-		"-- pulmtln --\n"
-		"per unit length capacitance and inductance for MTL cross - sections.\n"
-		"Visit https://github.com/OpenSEMBA/pulmtln for more information.\n"
+		"-- tulip --\n"
+		"Per unit length capacitance and inductance for MTL cross - sections.\n"
+		"Visit https://github.com/OpenSEMBA/tulip for more information.\n"
 		"Available options"
 	);
 	desc.add_options()
@@ -46,5 +46,5 @@ int main(int argc, char* argv[])
 
 	driver.run();
 
-	std::cout << "-- pulmtln finished succesfully --" << std::endl;
+	std::cout << "-- tulip finished succesfully --" << std::endl;
 }

src/pulmtln.cpp

@@ -221,7 +221,7 @@ TEST_F(DriverTest, five_wires)
 	ASSERT_EQ(couplingExpected.NumCols(), out.NumCols());
 
 	double rTol{ 0.05 };
-	for (int i{ 0 }; i < couplingExpected.NumRows(); i++) {
+for (int i{ 0 }; i < couplingExpected.NumRows(); i++) {
 		for (int j{ 0 }; j < couplingExpected.NumCols(); j++) {
 			EXPECT_LE(std::abs(couplingExpected(i, j) - out(i, j)), rTol)
 				<< "In C(" << i << ", " << j << ")";
@@ -725,6 +725,9 @@ TEST_F(DriverTest, lansink2024_single_wire_multipolar_in_cell_parameters)
 	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);
+
+	saveToJSONFile(inCell.toJSON(),
+		"lansink2024_single_wire_multipolar.inCellPotentials.out.json");
 }
 
 TEST_F(DriverTest, getCFromGeneralizedC_two_wires_open)
@@ -789,3 +792,141 @@ TEST_F(DriverTest, getCFromGeneralizedC_three_wires)
 		}
 	}
 }
+
+TEST_F(DriverTest, lansink2024_small_one_centered_fdtd_cell_vs_multipolar)
+{
+	// In-cell capacitances centered in conductor 0
+	// Using meshed FDTD cell.
+	InCellPotentials fdtdCellPotentials;
+	{
+		const std::string CASE{ "lansink2024_small_one_centered_fdtd_cell" };
+		fdtdCellPotentials = Driver::loadFromFile(
+			casesFolder() + CASE + "/" + CASE + ".pulmtln.in.json").getInCellPotentials();
+	}
+	auto fdtdCellComputedC00 = fdtdCellPotentials.getCapacitanceUsingInnerRegion(0, 0);
+	auto fdtdCellComputedC01 = fdtdCellPotentials.getCapacitanceUsingInnerRegion(0, 1);
+
+	// Using multipolar expansion.
+	InCellPotentials multipolarPotentials;
+	{
+		const std::string CASE{ "lansink2024_small_one_centered" };
+		multipolarPotentials = Driver::loadFromFile(
+			casesFolder() + CASE + "/" + CASE + ".pulmtln.in.json").getInCellPotentials();
+	}
+	Box fdtdCell{ {-0.1, -0.1}, {0.1, 0.1} };
+	auto multipolarComputedC00 = multipolarPotentials.getCapacitanceOnBox(0, 0, fdtdCell);
+	auto multipolarComputedC01 = multipolarPotentials.getCapacitanceOnBox(0, 1, fdtdCell);
+
+	// Compares results
+	EXPECT_NEAR(fdtdCellComputedC00, multipolarComputedC00, 0.03e-12);
+	EXPECT_NEAR(fdtdCellComputedC01, multipolarComputedC01, 0.03e-12);
+	EXPECT_LE(relError(fdtdCellComputedC00, multipolarComputedC00), 0.002);
+	EXPECT_LE(relError(fdtdCellComputedC01, multipolarComputedC01), 0.002);
+
+	saveToJSONFile(multipolarPotentials.toJSON(), 
+		"lansink2024_small_one_centered.inCellPotentials.out.json");
+}
+
+TEST_F(DriverTest, lansink2024_small_one_centered_different_integration_centers)
+{
+	// In-cell capacitances centered in conductor 0
+	InCellPotentials multipolarPotentials;
+		const std::string CASE{ "lansink2024_small_one_centered" };
+	
+	multipolarPotentials = Driver::loadFromFile(
+		casesFolder() + CASE + "/" + CASE + ".pulmtln.in.json").getInCellPotentials();
+
+	mfem::DenseMatrix geometricC(2,2);
+	Box fdtdCellCenteredOnConductor0{ {-0.1, -0.1}, {0.1, 0.1} };
+	{
+		geometricC(0, 0) = multipolarPotentials.getCapacitanceOnBox(0, 0, fdtdCellCenteredOnConductor0);
+		geometricC(0, 1) = multipolarPotentials.getCapacitanceOnBox(0, 1, fdtdCellCenteredOnConductor0);
+	}
+	{
+		Box fdtdCellCenteredOnConductor1{ {-0.12, -0.1}, {0.08, 0.1} };
+		geometricC(1, 0) = multipolarPotentials.getCapacitanceOnBox(1, 0, fdtdCellCenteredOnConductor1);
+		geometricC(1, 1) = multipolarPotentials.getCapacitanceOnBox(1, 1, fdtdCellCenteredOnConductor1);
+	}
+
+	mfem::DenseMatrix chargeCenteredC(2, 2);
+	for (int i = 0; i < 2; i++) {
+		Box chargeCenteredCell{ fdtdCellCenteredOnConductor0 };
+		chargeCenteredCell.displace(multipolarPotentials.electric.at(i).expansionCenter);
+		for (int j = 0; j < 2; j++) {
+			chargeCenteredC(i,j) = multipolarPotentials.getCapacitanceOnBox(i, j, chargeCenteredCell);
+		}
+	}
+	
+	// Compares results
+	for (int i = 0; i < 2; i++) {
+		for (int j = 0; j < 2; j++) {
+			EXPECT_NEAR(geometricC(i,j), chargeCenteredC(i,j), 0.1e-11);
+		}
+	}
+
+}
+
+TEST_F(DriverTest, DISABLED_realistic_case_with_dielectrics_fdtd_cell)
+{
+	InCellPotentials fdtdCellPotentials;
+	{
+		const std::string CASE{ "realistic_case_with_dielectrics_fdtd_cell" };
+		fdtdCellPotentials = Driver::loadFromFile(
+			casesFolder() + CASE + "/" + CASE + ".pulmtln.in.json").getInCellPotentials();
+	}
+	auto fdtdCellComputedC_0 = fdtdCellPotentials.getCapacitanceUsingInnerRegion(0, 0);
+	auto fdtdCellComputedC_16 = fdtdCellPotentials.getCapacitanceUsingInnerRegion(0, 16);
+	auto fdtdCellComputedC_25 = fdtdCellPotentials.getCapacitanceUsingInnerRegion(0, 25);
+	auto fdtdCellComputedC_30 = fdtdCellPotentials.getCapacitanceUsingInnerRegion(0, 30);
+
+	auto fdtdCellComputedL_0  = fdtdCellPotentials.getInductanceUsingInnerRegion(0, 0);
+	auto fdtdCellComputedL_16 = fdtdCellPotentials.getInductanceUsingInnerRegion(0, 16);
+	auto fdtdCellComputedL_25 = fdtdCellPotentials.getInductanceUsingInnerRegion(0, 25);
+	auto fdtdCellComputedL_30 = fdtdCellPotentials.getInductanceUsingInnerRegion(0, 30);
+}
+
+TEST_F(DriverTest, realistic_case_with_dielectrics)
+{
+	InCellPotentials mP;
+	{
+		const std::string CASE{ "realistic_case_with_dielectrics" };
+		mP = Driver::loadFromFile(
+			casesFolder() + CASE + "/" + CASE + ".pulmtln.in.json").getInCellPotentials();
+	}
+
+	Box fdtdCellCenteredOnConductor0{ {-0.016209, -0.009066}, {0.013791, 0.020934} };
+	auto mPComputedC_0  = mP.getCapacitanceOnBox(0,  0, fdtdCellCenteredOnConductor0);
+	auto mPComputedC_16 = mP.getCapacitanceOnBox(0, 16, fdtdCellCenteredOnConductor0);
+	auto mPComputedC_25 = mP.getCapacitanceOnBox(0, 25, fdtdCellCenteredOnConductor0);
+	auto mPComputedC_30 = mP.getCapacitanceOnBox(0, 30, fdtdCellCenteredOnConductor0);
+
+	auto mPComputedL_0  = mP.getInductanceOnBox(0,  0, fdtdCellCenteredOnConductor0);
+	auto mPComputedL_16 = mP.getInductanceOnBox(0, 16, fdtdCellCenteredOnConductor0);
+	auto mPComputedL_25 = mP.getInductanceOnBox(0, 25, fdtdCellCenteredOnConductor0);
+	auto mPComputedL_30 = mP.getInductanceOnBox(0, 30, fdtdCellCenteredOnConductor0);
+
+	// Compare with C and L computed using the fdtd cell.
+	const double rTol = 0.015;
+	const double fdtdCellComputed_C_0  = 4.0911726228481947e-11;
+	const double fdtdCellComputed_C_16 = 1.2547925523607968e-10;
+	const double fdtdCellComputed_C_25 = 5.9060595987059621e-11;
+	const double fdtdCellComputed_C_30 = 1.7797154919313720e-10;
+	const double fdtdCellComputed_L_0  = 2.9989786293920517e-07;
+	const double fdtdCellComputed_L_16 = 8.7458344767957649e-08;
+	const double fdtdCellComputed_L_25 = 2.0233814651758583e-07;
+	const double fdtdCellComputed_L_30 = 5.9647510363871327e-08;
+
+	EXPECT_LE(relError(fdtdCellComputed_C_0,  mPComputedC_0),  rTol) << "C(0,0) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_C_16, mPComputedC_16), rTol) << "C(0,16) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_C_25, mPComputedC_25), rTol) << "C(0,25) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_C_30, mPComputedC_30), rTol) << "C(0,30) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_L_0,  mPComputedL_0),  rTol) << "L(0,0) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_L_16, mPComputedL_16), rTol) << "L(0,16) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_L_25, mPComputedL_25), rTol) << "L(0,25) mismatch";
+	EXPECT_LE(relError(fdtdCellComputed_L_30, mPComputedL_30), rTol) << "L(0,30) mismatch";
+
+	 
+	// For debugging
+	saveToJSONFile(mP.toJSON(),
+		"realistic_case_with_dielectrics.inCellPotentials.out.json");
+}

test/DriverTest.cpp

@@ -0,0 +1,24 @@
+{
+   "geometries": [
+      {
+         "geometry": "Vacuum_0",
+         "area": 61787.60197630928
+      },
+      {
+         "geometry": "Dielectric_0",
+         "area": 39682.69914198743
+      },
+      {
+         "geometry": "Conductor_0",
+         "area": 3.141592653589793
+      },
+      {
+         "geometry": "Conductor_1",
+         "area": 314.1592653589793
+      },
+      {
+         "geometry": "OpenBoundary_0",
+         "area": 101787.60197630929
+      }
+   ]
+}