diff --git a/.vscode/launch.json b/.vscode/launch.json index 62b7839..2341c35 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -33,41 +33,11 @@ ] }, { - "name": "pulmtln_tests (gdb)", + "name": "tulip (gdb)", "type": "cppdbg", "request": "launch", - "program": "${workspaceFolder}/build-dbg/bin/pulmtln_tests", - "args": [], - "stopAtEntry": false, - "cwd": "${workspaceFolder}", - "environment": [], - "externalConsole": false, - "MIMode": "gdb", - "visualizerFile": [ - "${workspaceFolder}/resources/Eigen.natvis", - "${workspaceFolder}/resources/mfem.natvis", - "${workspaceFolder}/resources/nlohmann_json.natvis" - ], - "showDisplayString": true, - "setupCommands": [ - { - "description": "Enable pretty-printing for gdb", - "text": "-enable-pretty-printing", - "ignoreFailures": true - }, - { - "description": "Set Disassembly Flavor to Intel", - "text": "-gdb-set disassembly-flavor intel", - "ignoreFailures": true - } - ] - }, - { - "name": "pulmtln (gdb)", - "type": "cppdbg", - "request": "launch", - "program": "${workspaceFolder}/build-dbg/bin/pulmtln", - "args": [], + "program": "${workspaceFolder}/build-dbg/bin/tulip", + "args": ["-i", "${workspaceFolder}/testData/realistic_case_with_dielectrics/realistic_case_with_dielectrics_fdtd_cell.tulip.input.json"], "stopAtEntry": false, "cwd": "${workspaceFolder}", "environment": [], diff --git a/CMakePresets.json b/CMakePresets.json index d3928e6..f3405a6 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -44,13 +44,15 @@ "cacheVariables": { "TULIP_USE_MFEM_AS_SUBDIRECTORY": true, "MFEM_ENABLE_TESTING": false, - "MFEM_USE_OPENMP": false + "MFEM_USE_OPENMP": true, + "CMAKE_BUILD_TYPE": "Release" } }, { "name": "gnu-dbg", "displayName": "GNU g++ compiler", "inherits": "default", + "binaryDir": "${sourceDir}/build-dbg", "cacheVariables": { "TULIP_USE_MFEM_AS_SUBDIRECTORY": true, "MFEM_ENABLE_TESTING": false, diff --git a/scripts/compute_incell_parameters.py b/scripts/compute_incell_parameters.py index 13a0ac3..75e12c1 100644 --- a/scripts/compute_incell_parameters.py +++ b/scripts/compute_incell_parameters.py @@ -8,6 +8,13 @@ L[i,j] (H/m) for a fixed reference conductor i and all available conductors j in that material association. +The user must select one mode with --mode: + + * inner-region: uses getInCellCapacitanceUsingInnerRegion / + getInCellInductanceUsingInnerRegion formulas. + * on-box: uses getInCellCapacitanceOnBox / getInCellInductanceOnBox via + multipolar integration over the specified cell box. + Formulas (from InCellPotentials::getInCellCapacitanceUsingInnerRegion and getInCellInductanceUsingInnerRegion in Results.cpp): @@ -25,7 +32,8 @@ A_i = magnetic[j].conductorPotentials[i] Usage: - python scripts/compute_incell_parameters.py + python scripts/compute_incell_parameters.py --mode inner-region + python scripts/compute_incell_parameters.py --mode on-box --cell-box xmin ymin xmax ymax python scripts/compute_incell_parameters.py --run-tests """ @@ -41,6 +49,7 @@ # Physical constants (match tulip/src/driver/constants.h) EPSILON0_SI = 8.8541878176e-12 # F/m MU0_SI = 4.0e-7 * math.pi # H/m +GRID_INTEGRATION_SAMPLING_POINTS = 100 def load_json(json_path: str) -> dict: @@ -68,6 +77,140 @@ def inductance(magnetic_solutions: list, i: int, j: int) -> float: return (a_i - avg_a_j) / i_j * MU0_SI +def rel_error(expected: float, computed: float) -> float: + """Relative error used by DriverTests: |a-b| / max(|a|, |b|).""" + denom = max(abs(expected), abs(computed)) + if denom == 0.0: + return 0.0 + return abs(expected - computed) / denom + + +def make_box(xmin: float, ymin: float, xmax: float, ymax: float) -> dict: + """Build a Box-like dict using the C++ field names.""" + return {"min": [xmin, ymin], "max": [xmax, ymax]} + + +def box_area(box: dict) -> float: + """Compute area of a Box-like dict.""" + return (box["max"][0] - box["min"][0]) * (box["max"][1] - box["min"][1]) + + +def point_in_box(point: tuple, box: dict) -> bool: + """Check point inclusion using the same closed intervals as C++.""" + return ( + box["min"][0] <= point[0] <= box["max"][0] + and box["min"][1] <= point[1] <= box["max"][1] + ) + + +def build_integration_planes_for_box(integration_box: dict, inner_region_box: dict): + """Port of buildIntegrationPlanesForBox in Results.cpp.""" + for x in (0, 1): + if ( + integration_box["min"][x] > inner_region_box["min"][x] + or integration_box["max"][x] < inner_region_box["max"][x] + ): + raise ValueError("Integration box has to be larger than the inner region.") + + planes = [] + for x in (0, 1): + control_points = { + integration_box["min"][x], + integration_box["max"][x], + inner_region_box["min"][x], + inner_region_box["max"][x], + } + sorted_points = sorted(control_points) + + axis_planes = set() + for idx in range(1, len(sorted_points)): + prev = sorted_points[idx - 1] + curr = sorted_points[idx] + step = (curr - prev) / GRID_INTEGRATION_SAMPLING_POINTS + for k in range(GRID_INTEGRATION_SAMPLING_POINTS): + axis_planes.add(prev + k * step) + axis_planes.add(curr) + + planes.append(sorted(axis_planes)) + + return planes + + +def multipolar_expansion(position: tuple, ab: list, expansion_center: list) -> float: + """2D multipolar expansion matching multipolarExpansion.h.""" + rx = position[0] - expansion_center[0] + ry = position[1] - expansion_center[1] + r = math.hypot(rx, ry) + phi = math.atan2(ry, rx) + + res = 0.0 + for n, coeff in enumerate(ab): + an = coeff[0] + bn = coeff[1] + if n == 0: + res -= an * math.log(r) + else: + res += (an * math.cos(n * phi) + bn * math.sin(n * phi)) / (r ** n) + return res / (2.0 * math.pi) + + +def get_average_potential(solution: dict, inner_box: dict, outer_box: dict) -> float: + """Port of getAveragePotential in Results.cpp.""" + integration_planes = build_integration_planes_for_box(outer_box, inner_box) + outer_v = 0.0 + + for m in range(1, len(integration_planes[0])): + for n in range(1, len(integration_planes[1])): + x_min = integration_planes[0][m - 1] + x_max = integration_planes[0][m] + y_min = integration_planes[1][n - 1] + y_max = integration_planes[1][n] + + mid_point = (0.5 * (x_min + x_max), 0.5 * (y_min + y_max)) + area = (x_max - x_min) * (y_max - y_min) + if point_in_box(mid_point, inner_box): + continue + + outer_v += area * multipolar_expansion( + mid_point, + solution["ab"], + solution["expansionCenter"], + ) + + inner_v = solution["innerRegionAveragePotential"] * box_area(inner_box) + return (inner_v + outer_v) / box_area(outer_box) + + +def capacitance_on_box( + electric_solutions: list, + i: int, + j: int, + inner_region_box: dict, + cell_box: dict, +) -> float: + """Compute C[i,j] in F/m using getInCellCapacitanceOnBox semantics.""" + sol_j = electric_solutions[j] + q_j = sol_j["ab"][0][0] + avg_v_j = get_average_potential(sol_j, inner_region_box, cell_box) + v_i = sol_j["conductorPotentials"][i] + return q_j / (v_i - avg_v_j) * EPSILON0_SI + + +def inductance_on_box( + magnetic_solutions: list, + i: int, + j: int, + inner_region_box: dict, + cell_box: dict, +) -> float: + """Compute L[i,j] in H/m using getInCellInductanceOnBox semantics.""" + sol_j = magnetic_solutions[j] + i_j = sol_j["ab"][0][0] + avg_a_j = get_average_potential(sol_j, inner_region_box, cell_box) + a_i = sol_j["conductorPotentials"][i] + return (a_i - avg_a_j) / i_j * MU0_SI + + def local_index_from_element_id(element_ids: list, element_id: int): """Map a physical conductor id to the local multipolar solution index.""" try: @@ -84,17 +227,26 @@ def element_ids_for_material(material_id: int, mat_assoc: list, fallback_size: i return list(range(fallback_size)) -def compute_material_rows(material: dict, mat_assoc: list, i_ref: int): - """Compute (j, C[i_ref,j], L[i_ref,j]) rows for one material.""" +def compute_material_rows(material: dict, mat_assoc: list, i_ref: int, mode: str, cell_box: dict = None): + """Compute (j, C[i_ref,j], L[i_ref,j]) rows for one material and mode.""" mat_id = material.get("id") mat_type = material.get("type", "unknown") mp = material.get("multipolarExpansion", {}) e_sols = mp.get("electric", []) m_sols = mp.get("magnetic", []) + inner_region_box = mp.get("innerRegionBox") if not e_sols or not m_sols: return None, "missing electric or magnetic solutions" + if mode == "on-box": + if cell_box is None: + return None, "on-box mode requires a cell box" + if inner_region_box is None: + return None, "on-box mode requires multipolarExpansion.innerRegionBox" + if "expansionCenter" not in e_sols[0] or "expansionCenter" not in m_sols[0]: + return None, "on-box mode requires expansionCenter in electric/magnetic solutions" + element_ids = element_ids_for_material(mat_id, mat_assoc, len(e_sols)) i_local = local_index_from_element_id(element_ids, i_ref) if i_local is None: @@ -105,8 +257,14 @@ def compute_material_rows(material: dict, mat_assoc: list, i_ref: int): j_local = local_index_from_element_id(element_ids, j) if j_local is None: continue - c_val = capacitance(e_sols, i_local, j_local) - l_val = inductance(m_sols, i_local, j_local) + if mode == "inner-region": + c_val = capacitance(e_sols, i_local, j_local) + l_val = inductance(m_sols, i_local, j_local) + elif mode == "on-box": + c_val = capacitance_on_box(e_sols, i_local, j_local, inner_region_box, cell_box) + l_val = inductance_on_box(m_sols, i_local, j_local, inner_region_box, cell_box) + else: + return None, f"unsupported mode '{mode}'" rows.append((j, c_val, l_val)) return { @@ -114,6 +272,7 @@ def compute_material_rows(material: dict, mat_assoc: list, i_ref: int): "mat_type": mat_type, "element_ids": element_ids, "rows": rows, + "mode": mode, }, None @@ -122,6 +281,7 @@ def print_material_report(report: dict, i_ref: int): print(f"\n{'=' * 55}") print(f" Material id={report['mat_id']} type={report['mat_type']}") print(f" Total conductors: {len(report['element_ids'])}") + print(f" Mode: {report['mode']}") print(f"{'=' * 55}") print(f" {'j':>4} {f'C[{i_ref},j] (F/m)':>18} {f'L[{i_ref},j] (H/m)':>18}") print(f" {'-' * 4} {'-' * 18} {'-' * 18}") @@ -140,6 +300,18 @@ def parse_args(argv: list): action="store_true", help="Run the unit tests embedded in this file.", ) + parser.add_argument( + "--mode", + choices=["inner-region", "on-box"], + help="Choose the computational mode. Required unless --run-tests is used.", + ) + parser.add_argument( + "--cell-box", + nargs=4, + type=float, + metavar=("XMIN", "YMIN", "XMAX", "YMAX"), + help="Cell box for on-box mode: xmin ymin xmax ymax.", + ) parser.add_argument( "json_path", nargs="?", @@ -156,8 +328,20 @@ def parse_args(argv: list): if args.run_tests: return args + if args.mode is None: + parser.error("--mode is required unless --run-tests is used") if args.json_path is None or args.i_ref is None: parser.error("json_path and i_ref are required unless --run-tests is used") + + if args.mode == "on-box": + if args.cell_box is None: + parser.error("--cell-box is required when --mode on-box is selected") + if args.cell_box[2] <= args.cell_box[0] or args.cell_box[3] <= args.cell_box[1]: + parser.error("--cell-box must satisfy xmax > xmin and ymax > ymin") + else: + if args.cell_box is not None: + parser.error("--cell-box is only valid when --mode on-box is selected") + return args @@ -175,9 +359,16 @@ def run(argv=None) -> int: json_path = os.path.abspath(args.json_path) i_ref = args.i_ref + mode = args.mode + cell_box = None + if mode == "on-box": + cell_box = make_box(*args.cell_box) print(f"Output JSON path : {json_path}") print(f"Reference conductor i : {i_ref}") + print(f"Mode : {mode}") + if cell_box is not None: + print(f"Cell box : {cell_box}") data = load_json(json_path) materials = data.get("materials", []) @@ -188,7 +379,7 @@ def run(argv=None) -> int: return 0 for material in materials: - report, warn = compute_material_rows(material, mat_assoc, i_ref) + report, warn = compute_material_rows(material, mat_assoc, i_ref, mode, cell_box) if warn is not None: print(f"[warn] Material {material.get('id')}: {warn}, skipping.") continue @@ -239,7 +430,7 @@ def test_compute_material_rows_uses_all_conductors(self): } mat_assoc = [{"materialId": 7, "elementIds": [16, 30]}] - report, warn = compute_material_rows(material, mat_assoc, 16) + report, warn = compute_material_rows(material, mat_assoc, 16, mode="inner-region") self.assertIsNone(warn) self.assertIsNotNone(report) self.assertEqual([row[0] for row in report["rows"]], [16, 30]) @@ -301,7 +492,7 @@ def test_reference_case_values_from_embedded_fixture(self): } mat_assoc = [{"elementIds": [16, 30], "materialId": 1}] - report, warn = compute_material_rows(material, mat_assoc, 16) + report, warn = compute_material_rows(material, mat_assoc, 16, mode="inner-region") self.assertIsNone(warn) self.assertIsNotNone(report) @@ -314,15 +505,104 @@ def test_reference_case_values_from_embedded_fixture(self): self.assertAlmostEqual(rows[30][0], 5.165814539740492e-11, delta=1e-20) self.assertAlmostEqual(rows[30][1], 2.1538714707844526e-07, delta=1e-15) - def test_parse_args_requires_positionals(self): + def test_on_box_mode_matches_driver_test_single_wire_values(self): + # This fixture is equivalent to the in-cell output used by + # DriverTest::lansink2024_single_wire_multipolar_in_cell_parameters. + material = { + "id": 1, + "type": "unshieldedMultiwire", + "multipolarExpansion": { + "innerRegionBox": { + "max": [0.0040000002, 0.0040000002], + "min": [-0.0040000002, -0.0040000002], + }, + "electric": [ + { + "ab": [ + [0.9766734089848975, 0.0], + [-5.686345063943916e-12, -3.090752789506213e-12], + [4.637444003229173e-11, 8.593758249532865e-12], + [7.014114840799095e-15, 4.218755938827563e-15], + [-9.241748384573654e-18, 3.9674592343454175e-18], + [-1.0388931187475782e-21, 3.5460599095173054e-21], + ], + "conductorPotentials": [1.0], + "expansionCenter": [-4.899147920455392e-08, -2.6628800768531107e-08], + "innerRegionAveragePotential": 0.9040784423949009, + } + ], + "magnetic": [ + { + "ab": [ + [0.9092956935239767, 0.0], + [4.676064105592532e-12, -1.9984633243172434e-12], + [1.79380465488009e-11, 3.6833862710835886e-12], + [-2.3149425432828747e-15, 1.8113967798566277e-15], + [8.898740430780706e-19, -6.653588044528165e-19], + [-4.9778872535021934e-23, -2.6188498558876825e-22], + ], + "conductorPotentials": [1.0], + "expansionCenter": [4.3634340291863276e-08, -1.8648510194811956e-08], + "innerRegionAveragePotential": 0.8490379205671101, + } + ], + }, + } + mat_assoc = [{"materialId": 1, "elementIds": [0]}] + fdtd_cell = make_box(-0.0075, -0.0075, 0.0075, 0.0075) + + report, warn = compute_material_rows( + material, + mat_assoc, + 0, + mode="on-box", + cell_box=fdtd_cell, + ) + self.assertIsNone(warn) + self.assertIsNotNone(report) + + rows = {j: (c_val, l_val) for j, c_val, l_val in report["rows"]} + self.assertIn(0, rows) + + expected_c00 = 49.11e-12 + expected_l00 = 320e-9 + r_tol = 0.06 + + self.assertLess(rel_error(expected_c00, rows[0][0]), r_tol) + self.assertLess(rel_error(expected_l00, rows[0][1]), r_tol) + + def test_parse_args_requires_mode(self): with contextlib.redirect_stderr(io.StringIO()): with self.assertRaises(SystemExit): - parse_args([]) + parse_args(["case.tulip.out.json", "16"]) - def test_parse_args_happy_path(self): - args = parse_args(["case.tulip.out.json", "16"]) + def test_parse_args_inner_region_happy_path(self): + args = parse_args(["--mode", "inner-region", "case.tulip.out.json", "16"]) self.assertEqual(args.json_path, "case.tulip.out.json") self.assertEqual(args.i_ref, 16) + self.assertEqual(args.mode, "inner-region") + + def test_parse_args_on_box_requires_cell_box(self): + with contextlib.redirect_stderr(io.StringIO()): + with self.assertRaises(SystemExit): + parse_args(["--mode", "on-box", "case.tulip.out.json", "16"]) + + def test_parse_args_on_box_happy_path(self): + args = parse_args( + [ + "--mode", + "on-box", + "--cell-box", + "-0.1", + "-0.1", + "0.1", + "0.1", + "case.tulip.out.json", + "16", + ] + ) + self.assertEqual(args.mode, "on-box") + self.assertEqual(args.cell_box, [-0.1, -0.1, 0.1, 0.1]) if __name__ == "__main__": diff --git a/src/Launcher.cpp b/src/Launcher.cpp index 72c230a..9fd83c8 100644 --- a/src/Launcher.cpp +++ b/src/Launcher.cpp @@ -4,8 +4,10 @@ #include "AdaptedInputParser.h" #include "Driver.h" +#include #include #include +#include #include #include @@ -13,6 +15,13 @@ namespace tulip { namespace { +using Clock = std::chrono::steady_clock; + +double elapsedSeconds(Clock::time_point start) +{ + return std::chrono::duration(Clock::now() - start).count(); +} + bool hasSuffix(const std::string& value, const std::string& suffix) { return value.size() >= suffix.size() && @@ -49,6 +58,28 @@ std::string extractCaseName(const std::string& inputFile) return std::filesystem::path(filename).stem().string(); } +void printTimingSummary( + bool hasMeshTiming, + double meshSeconds, + const DriverTimings& timings) +{ + std::cout << "Timing summary:" << std::endl; + std::cout << std::fixed << std::setprecision(3); + if (hasMeshTiming) { + std::cout << " mesh: " << meshSeconds << " s" << std::endl; + } + else { + std::cout << " mesh: N/A" << std::endl; + } + std::cout << " solve: " << timings.solveSeconds << " s" << std::endl; + if (timings.multipolarComputed) { + std::cout << " multipolar: " << timings.multipolarSeconds << " s" << std::endl; + } + else { + std::cout << " multipolar: N/A" << std::endl; + } +} + } // namespace Launcher::Launcher(const std::string& inputFile, const std::string& exportFolder) @@ -65,12 +96,15 @@ void Launcher::run() std::cout << "Loading input file: " << inputFile_ << std::endl; const std::string caseNamePrefix = extractCaseName(inputFile_) + "."; const std::string outputPathPrefix = ensureTrailingSlash(exportFolder_) + caseNamePrefix; + bool hasMeshTiming = false; + double meshSeconds = 0.0; if (isAdaptedJson(inputFile_)) { auto driver = Driver::loadFromAdaptedFile(inputFile_); driver.setExportFolder(outputPathPrefix); std::cout << "Running Tulip analysis..." << std::endl; driver.run(); + printTimingSummary(hasMeshTiming, meshSeconds, driver.getTimings()); } else if (isInputJson(inputFile_)) { const bool initializedHere = !gmsh::isInitialized(); @@ -79,12 +113,16 @@ void Launcher::run() } try { + const auto meshStart = Clock::now(); Adapter adapter(inputFile_); + meshSeconds = elapsedSeconds(meshStart); + hasMeshTiming = true; AdaptedInputParser parser(inputFile_, adapter.getAdaptedInputJSON()); Driver driver(parser.readModel(), parser.readDriverOptions()); driver.setExportFolder(outputPathPrefix); std::cout << "Running Tulip analysis..." << std::endl; driver.run(); + printTimingSummary(hasMeshTiming, meshSeconds, driver.getTimings()); } catch (...) { if (initializedHere) { diff --git a/src/adapter/Adapter.cpp b/src/adapter/Adapter.cpp index b57108c..7dd42db 100644 --- a/src/adapter/Adapter.cpp +++ b/src/adapter/Adapter.cpp @@ -16,15 +16,6 @@ namespace tulip { -const Adapter::MeshingOptions Adapter::DEFAULT_MESHING_OPTIONS{ - {"Mesh.MshFileVersion", 2.2}, - {"Mesh.MeshSizeFromCurvature", 50.0}, - {"Mesh.ElementOrder", 3.0}, - {"Mesh.ScalingFactor", 1e-3}, - {"Mesh.SurfaceFaces", 1.0}, - {"Mesh.MeshSizeMax", 40.0}, -}; - std::pair, std::vector>> @@ -799,7 +790,6 @@ nlohmann::json buildAdaptedJson(const std::string& caseName, } // namespace Adapter::Adapter(const std::string& inputFile) - : allShapes(EntityList{}, nlohmann::json::object()) { const std::filesystem::path inputPath(inputFile); if (!hasSuffix(inputPath.filename().string(), ".tulip.input.json")) { @@ -813,7 +803,6 @@ Adapter::Adapter(const std::string& inputFile) Adapter::Adapter(const nlohmann::json& inputJson, const std::string& caseName, const std::string& inputDir) - : allShapes(EntityList{}, inputJson) { initialize(inputJson, caseName, inputDir); } @@ -846,7 +835,11 @@ void Adapter::initialize(const nlohmann::json& inputJson, gmsh::model::occ::synchronize(); validateLayerNamesMatchStep(inputJson, shapes); - allShapes = ShapesClassification(shapes, inputJson); + allShapes = ShapesClassification( + shapes, + inputJson, + adapterOptions_.innerRegionBoxScalingFactor, + adapterOptions_.farRegionDiskScalingFactor); allShapes.ensureDielectricsDoNotOverlap(); allShapes.vacuum = allShapes.buildVacuumDomain(); @@ -874,7 +867,6 @@ void Adapter::initialize(const nlohmann::json& inputJson, for (const auto& [opt, val] : adapterOptions_.gmshOptions) { gmsh::option::setNumber(opt, val); } - gmsh::model::mesh::generate(2); gmsh::model::mesh::removeDuplicateNodes(); diff --git a/src/adapter/Adapter.h b/src/adapter/Adapter.h index f7fa7f8..fa28032 100644 --- a/src/adapter/Adapter.h +++ b/src/adapter/Adapter.h @@ -17,9 +17,6 @@ class Adapter { public: using MeshingOptions = std::map; - // Baseline gmsh options used by adapter tests and CLI code. - static const MeshingOptions DEFAULT_MESHING_OPTIONS; - Adapter(const std::string& inputFile); Adapter(const nlohmann::json& inputJson, const std::string& caseName, diff --git a/src/adapter/AdapterOptions.h b/src/adapter/AdapterOptions.h index 8f0a673..3da076b 100644 --- a/src/adapter/AdapterOptions.h +++ b/src/adapter/AdapterOptions.h @@ -10,7 +10,7 @@ struct AdapterOptions // Multiplication factor applied to the size of the box created // to determine the inner region used for unshielded multiwires // bundles with no open boundary defined. - double innerRegionBoxScalingFactor = 1.15; + double innerRegionBoxScalingFactor = 1.30; // For an open case in which no open boundary is defined a // circle is created enclosing all the defined entities. @@ -25,11 +25,11 @@ struct AdapterOptions // User specified options modify these options. std::map gmshOptions = { {"Mesh.MshFileVersion", 2.2}, - {"Mesh.MeshSizeFromCurvature", 50.0}, + {"Mesh.MeshSizeFromCurvature", 40.0}, {"Mesh.ElementOrder", 3.0}, {"Mesh.ScalingFactor", 1e-3}, {"Mesh.SurfaceFaces", 1.0}, - {"Mesh.MeshSizeMax", 40.0} + {"Mesh.MeshSizeMax", 2.5} }; }; diff --git a/src/adapter/ShapesClassification.cpp b/src/adapter/ShapesClassification.cpp index dd11ed6..719b513 100644 --- a/src/adapter/ShapesClassification.cpp +++ b/src/adapter/ShapesClassification.cpp @@ -11,8 +11,6 @@ namespace tulip { namespace { -constexpr double innerRegionBoxScalingFactor = 1.15; -constexpr double farRegionBoxScalingFactor = 4.0; constexpr double defaultRelativePermittivity = 1.0; constexpr double relativePermittivityTolerance = 1e-12; constexpr double intersectionAreaTolerance = 1e-18; @@ -164,7 +162,9 @@ bool anyPairOfConductorsOverlap(const EntityMap& conductors) } ShapesClassification::ShapesClassification(const EntityList& shapes, - const std::string& jsonFile) + const std::string& jsonFile, + double innerRegionBoxScalingFactor, + double farRegionDiskScalingFactor) : ShapesClassification(shapes, [&jsonFile]() { std::ifstream f(jsonFile); if (!f.is_open()) { @@ -173,11 +173,15 @@ ShapesClassification::ShapesClassification(const EntityList& shapes, nlohmann::json jsonData; f >> jsonData; return jsonData; - }()) + }(), innerRegionBoxScalingFactor, farRegionDiskScalingFactor) {} ShapesClassification::ShapesClassification(const EntityList& shapes, - const nlohmann::json& jsonData) + const nlohmann::json& jsonData, + double innerRegionBoxScalingFactor, + double farRegionDiskScalingFactor) + : innerRegionBoxScalingFactor_(innerRegionBoxScalingFactor), + farRegionDiskScalingFactor_(farRegionDiskScalingFactor) { gmsh::model::occ::synchronize(); @@ -475,7 +479,7 @@ EntityMap ShapesClassification::buildOpenVacuumDomain() { auto lengths = boundingBox.getLengths(); double bbMaxLen = *std::max_element(lengths.begin(), lengths.end()); - double nearBoxSize = bbMaxLen * innerRegionBoxScalingFactor; + double nearBoxSize = bbMaxLen * innerRegionBoxScalingFactor_; auto center = boundingBox.getCenter(); double nVX = center[0] - nearBoxSize / 2.0; double nVY = center[1] - nearBoxSize / 2.0; @@ -485,7 +489,7 @@ EntityMap ShapesClassification::buildOpenVacuumDomain() { nVX, nVY, nVZ, nearBoxSize, nearBoxSize); EntityList nearVacuum = {{2, nearRectTag}}; - double farDiameter = farRegionBoxScalingFactor * boundingBox.getDiagonal(); + double farDiameter = farRegionDiskScalingFactor_ * boundingBox.getDiagonal(); int farDiskTag = gmsh::model::occ::addDisk( center[0], center[1], center[2], farDiameter, farDiameter); EntityList farVacuum = {{2, farDiskTag}}; diff --git a/src/adapter/ShapesClassification.h b/src/adapter/ShapesClassification.h index 8048c4a..85d51d6 100644 --- a/src/adapter/ShapesClassification.h +++ b/src/adapter/ShapesClassification.h @@ -27,8 +27,15 @@ class ShapesClassification { EntityList allShapes; Graph nestedGraph; - ShapesClassification(const EntityList& shapes, const std::string& jsonFile); - ShapesClassification(const EntityList& shapes, const nlohmann::json& jsonData); + ShapesClassification() = default; + ShapesClassification(const EntityList& shapes, + const std::string& jsonFile, + double innerRegionBoxScalingFactor, + double farRegionDiskScalingFactor); + ShapesClassification(const EntityList& shapes, + const nlohmann::json& jsonData, + double innerRegionBoxScalingFactor, + double farRegionDiskScalingFactor); static int getNumberFromName(const std::string& entityName, const std::string& label); @@ -53,6 +60,8 @@ class ShapesClassification { private: nlohmann::json crossSectionData_; + double innerRegionBoxScalingFactor_; + double farRegionDiskScalingFactor_; double getDielectricRelativePermittivity(const std::string& geometryName) const; bool dielectricHasPriorityOver(const std::string& lhs, diff --git a/src/driver/AdaptedInputParser.cpp b/src/driver/AdaptedInputParser.cpp index 038c19b..2d43b87 100644 --- a/src/driver/AdaptedInputParser.cpp +++ b/src/driver/AdaptedInputParser.cpp @@ -166,6 +166,7 @@ DriverOptions Parser::readDriverOptions() const const auto& j = json_.at("driverOptions"); DriverOptions res; + setIfExists(j, res.multipolarExpansionOrder, "multipolarExpansionOrder"); setIfExists(j, res.solverOptions.order, "order"); setIfExists(j, res.solverOptions.printIterations, "printIterations"); diff --git a/src/driver/Driver.cpp b/src/driver/Driver.cpp index c7fca3b..eafa8d4 100644 --- a/src/driver/Driver.cpp +++ b/src/driver/Driver.cpp @@ -4,12 +4,21 @@ #include "AdaptedInputParser.h" #include "multipolarExpansion.h" +#include + using namespace mfem; namespace tulip { namespace { +using Clock = std::chrono::steady_clock; + +double elapsedSeconds(Clock::time_point start) +{ + return std::chrono::duration(Clock::now() - start).count(); +} + void configureMfemDeviceIfAvailable() { #ifdef MFEM_USE_OPENMP @@ -123,6 +132,8 @@ Driver::Driver(Model&& model, const DriverOptions& opts) : throw std::runtime_error("Model must have at least one conductor."); } + const auto solveStart = Clock::now(); + // Solve for all conductors. std::cout << "Solving electrostatic problems:" << std::endl; electric_ = solveForAllConductors(FieldType::electric); @@ -134,6 +145,8 @@ Driver::Driver(Model&& model, const DriverOptions& opts) : else { std::cout << "No dielectrics found. Reusing electrostatic solution for magnetostatic." << std::endl; } + + timings_.solveSeconds = elapsedSeconds(solveStart); } mfem::DenseMatrix Driver::getCFromGeneralizedC( @@ -364,10 +377,6 @@ MultiwireParametersByDomain Driver::getMultiwireParametersByDomains() DomainTree domainTree{ idToDomain }; res.setDomainTree(domainTree); const auto openness = model_.determineOpenness(); - std::unique_ptr openDomainPotentials; - if (openness == Model::Openness::open) { - openDomainPotentials = std::make_unique(getInCellPotentials()); - } // Build mappings from conductor ID to model data. auto conductors = model_.getMaterials().getConductors(); @@ -383,6 +392,12 @@ MultiwireParametersByDomain Driver::getMultiwireParametersByDomains() auto globalGC = getGeneralizedCMatrix(FieldType::electric); auto globalGC0 = getGeneralizedCMatrix(FieldType::magnetic); + std::unique_ptr openDomainPotentials; + if (openness == Model::Openness::open) { + openDomainPotentials = std::make_unique( + getInCellPotentials(&globalGC, &globalGC0)); + } + for (const auto& [domId, domain] : idToDomain) { if (openness == Model::Openness::open && domain.ground == Domain::UNDEFINED_GROUND) { @@ -655,7 +670,8 @@ void Driver::loadFloatingPotentials( } std::map Driver::getFieldParameters( - FieldType fieldType) + FieldType fieldType, + const mfem::DenseMatrix* generalizedC) { std::map res; @@ -667,16 +683,21 @@ std::map Driver::getFieldParameters( << " field coefficients." << std::endl; const auto conductors = model_.getMaterials().getConductors(); - // Compute the C matrix once for all conductors to avoid - // reassembling operators N times (was O(N^2), now O(N)). - mfem::DenseMatrix C = getGeneralizedCMatrix(fieldType); + std::unique_ptr computedC; + if (generalizedC == nullptr) { + // Compute the C matrix once for all conductors to avoid + // reassembling operators N times (was O(N^2), now O(N)). + computedC = std::make_unique( + getGeneralizedCMatrix(fieldType)); + generalizedC = computedC.get(); + } for (const auto& cI : conductors) { auto condI = cI->getConductorId(); std::cout << "- Conductor #" << condI << "... " << std::flush; - auto fp = computeFloatingPotentialsFromC(condI, C); + auto fp = computeFloatingPotentialsFromC(condI, *generalizedC); loadFloatingPotentials(sP, fp); @@ -688,7 +709,9 @@ std::map Driver::getFieldParameters( std::copy( centerOfCharge.begin(), centerOfCharge.end(), res[condI].expansionCenter.begin()); - res[condI].ab = s.getMultipolarCoefficients(opts_.multipolarExpansionOrder); + res[condI].ab = s.getMultipolarCoefficients( + opts_.multipolarExpansionOrder, + centerOfCharge); for (const auto& cJ : conductors) { auto condJ = cJ->getConductorId(); res[condI].conductorPotentials[condJ] = fp.at(condJ); @@ -700,18 +723,34 @@ std::map Driver::getFieldParameters( } InCellPotentials Driver::getInCellPotentials() +{ + return getInCellPotentials(nullptr, nullptr); +} + +InCellPotentials Driver::getInCellPotentials( + const mfem::DenseMatrix* electricGeneralizedC, + const mfem::DenseMatrix* magneticGeneralizedC) { if (model_.determineOpenness() != Model::Openness::open) { throw std::runtime_error("In cell potentials can only be determined for open problems."); } + const auto multipolarStart = Clock::now(); + InCellPotentials res; res.getInnerRegionBox() = model_.getInnerRegionBoundingBox(); - res.getElectric() = getFieldParameters(FieldType::electric); - res.getMagnetic() = getFieldParameters(FieldType::magnetic); + res.getElectric() = getFieldParameters( + FieldType::electric, + electricGeneralizedC); + res.getMagnetic() = getFieldParameters( + FieldType::magnetic, + magneticGeneralizedC); + + timings_.multipolarSeconds += elapsedSeconds(multipolarStart); + timings_.multipolarComputed = true; return res; } -} \ No newline at end of file +} diff --git a/src/driver/Driver.h b/src/driver/Driver.h index a3dc835..c29452a 100644 --- a/src/driver/Driver.h +++ b/src/driver/Driver.h @@ -14,6 +14,12 @@ struct SolvedProblem { std::map solutions; }; +struct DriverTimings { + double solveSeconds{ 0.0 }; + double multipolarSeconds{ 0.0 }; + bool multipolarComputed{ false }; +}; + class Driver { public: enum class FieldType { @@ -47,6 +53,7 @@ class Driver { DenseMatrix getLMatrix(); const Model& getModel() const { return model_; } + const DriverTimings& getTimings() const { return timings_; } static Driver loadFromAdaptedFile(const std::string& filename); static Driver adaptFromFile(const std::string& filename); @@ -63,6 +70,7 @@ class Driver { Model model_; DriverOptions opts_; SolvedProblem electric_, magnetic_; + DriverTimings timings_; SolvedProblem solveForAllConductors( FieldType fieldType); @@ -71,8 +79,12 @@ class Driver { double getInnerRegionAveragePotential( const Solver& s, bool includeConductors); + InCellPotentials getInCellPotentials( + const mfem::DenseMatrix* electricGeneralizedC, + const mfem::DenseMatrix* magneticGeneralizedC); std::map getFieldParameters( - FieldType fieldType); + FieldType fieldType, + const mfem::DenseMatrix* generalizedC = nullptr); std::map computeFloatingPotentialsFromC( ConductorId prescribedId, const mfem::DenseMatrix& C) const; @@ -80,4 +92,4 @@ class Driver { }; -} \ No newline at end of file +} diff --git a/src/driver/Solver.cpp b/src/driver/Solver.cpp index f6b29ee..2f416a4 100644 --- a/src/driver/Solver.cpp +++ b/src/driver/Solver.cpp @@ -372,13 +372,19 @@ multipolarCoefficients Solver::getMultipolarCoefficients( std::size_t order) const { auto centerOfCharge{ getCenterOfCharge() }; + return getMultipolarCoefficients(order, centerOfCharge); +} +multipolarCoefficients Solver::getMultipolarCoefficients( + std::size_t order, + const Vector& center) const +{ multipolarCoefficients ab(order + 1); for (int n = 0; n < order + 1; n++) { ab[n] = { - getChargeMomentComponent(n, 0, centerOfCharge), - getChargeMomentComponent(n, 1, centerOfCharge) + getChargeMomentComponent(n, 0, center), + getChargeMomentComponent(n, 1, center) }; } @@ -542,4 +548,4 @@ void Solver::writeParaViewFields(ParaViewDataCollection& pv) const pv.Save(); } -} \ No newline at end of file +} diff --git a/src/driver/Solver.h b/src/driver/Solver.h index 96b11a4..f4cf0c2 100644 --- a/src/driver/Solver.h +++ b/src/driver/Solver.h @@ -49,6 +49,9 @@ class Solver { Vector getCenterOfCharge() const; double getChargeMomentComponent(int order, int component, const Vector& center) const; multipolarCoefficients getMultipolarCoefficients(std::size_t order) const; + multipolarCoefficients getMultipolarCoefficients( + std::size_t order, + const Vector& center) const; SolverSolution getSolution() const; void setSolution(const SolverSolution&); @@ -99,4 +102,4 @@ class Solver { GridFunction& gf) const; }; -} \ No newline at end of file +} diff --git a/test/adapter/AdapterTest.cpp b/test/adapter/AdapterTest.cpp index 0fac84a..2a544de 100644 --- a/test/adapter/AdapterTest.cpp +++ b/test/adapter/AdapterTest.cpp @@ -13,6 +13,7 @@ #include "Mesher.h" #include "ShapesClassification.h" #include "TestUtils.h" +#include "AdapterOptions.h" using namespace tulip; @@ -410,7 +411,14 @@ TEST_F(AdapterTest, overlapping_dielectrics_prioritize_higher_relative_permittiv })} }; - ShapesClassification classification(shapes, inputJson); + AdapterOptions opts; + + ShapesClassification classification( + shapes, + inputJson, + opts.innerRegionBoxScalingFactor, + opts.farRegionDiskScalingFactor + ); classification.ensureDielectricsDoNotOverlap(); gmsh::model::occ::synchronize(); @@ -440,7 +448,14 @@ TEST_F(AdapterTest, shapes_classification_without_roots_is_treated_as_open_probl {"layers", nlohmann::json::array()} }; - ShapesClassification classification(shapes, inputJson); + AdapterOptions opts; + + ShapesClassification classification( + shapes, + inputJson, + opts.innerRegionBoxScalingFactor, + opts.farRegionDiskScalingFactor + ); EXPECT_TRUE(classification.isOpenCase); EXPECT_TRUE(classification.isOpenProblem()); diff --git a/test/driver/ParserTest.cpp b/test/driver/ParserTest.cpp index 6ea0699..e5fe12a 100644 --- a/test/driver/ParserTest.cpp +++ b/test/driver/ParserTest.cpp @@ -16,6 +16,7 @@ TEST_F(ParserTest, empty_coax) }; auto opts{ parser.readDriverOptions() }; + EXPECT_EQ(3, opts.multipolarExpansionOrder); EXPECT_EQ(3, opts.solverOptions.order); EXPECT_EQ(true, opts.exportParaViewSolution); @@ -52,6 +53,7 @@ TEST_F(ParserTest, partially_filled_coax) }; auto opts{ parser.readDriverOptions() }; + EXPECT_EQ(3, opts.multipolarExpansionOrder); EXPECT_EQ(3, opts.solverOptions.order); EXPECT_EQ(true, opts.exportParaViewSolution); @@ -106,3 +108,25 @@ TEST_F(ParserTest, hasDielectrics_returns_true_when_dielectric_present) EXPECT_TRUE(model.getMaterials().hasDielectrics()); } +TEST_F(ParserTest, reads_custom_multipolar_expansion_order) +{ + nlohmann::json adaptedJson = { + {"driverOptions", { + {"multipolarExpansionOrder", 5}, + {"order", 2}, + {"printIterations", true}, + {"exportParaViewSolution", false}, + {"exportFolder", "Results/custom/"} + }} + }; + + Parser parser{ "unused.tulip.adapted.json", adaptedJson }; + auto opts{ parser.readDriverOptions() }; + + EXPECT_EQ(5, opts.multipolarExpansionOrder); + EXPECT_EQ(2, opts.solverOptions.order); + EXPECT_EQ(true, opts.solverOptions.printIterations); + EXPECT_EQ(false, opts.exportParaViewSolution); + EXPECT_EQ("Results/custom/", opts.exportFolder); +} + diff --git a/testData/realistic_case_with_dielectrics/realistic_case_with_dielectrics_inner_region.tulip.input.json b/testData/realistic_case_with_dielectrics/realistic_case_with_dielectrics_inner_region.tulip.input.json index 5ba643f..9cf4cb5 100644 --- a/testData/realistic_case_with_dielectrics/realistic_case_with_dielectrics_inner_region.tulip.input.json +++ b/testData/realistic_case_with_dielectrics/realistic_case_with_dielectrics_inner_region.tulip.input.json @@ -1,7 +1,8 @@ { "driverOptions": { + "multipolarExpansionOrder": 3, "exportFolder": "Results/realistic_case_with_dielectrics_inner_region/", - "exportParaViewSolution": true + "exportParaViewSolution": false }, "materials": [ {"id": 1, "type": "conductor"},