From 33bcbe15285023c5bd4e80fc7ff7ffdf2fddee21 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 5 May 2026 11:30:44 +0200 Subject: [PATCH 1/3] Adds script to compute incell from outputs. --- scripts/compute_incell_parameters.py | 139 +++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 scripts/compute_incell_parameters.py diff --git a/scripts/compute_incell_parameters.py b/scripts/compute_incell_parameters.py new file mode 100644 index 0000000..1486ad1 --- /dev/null +++ b/scripts/compute_incell_parameters.py @@ -0,0 +1,139 @@ +""" +compute_incell_parameters.py +============================= +Postprocessing script for a Tulip solved case. + +Reads the .tulip.out.json file and, for each material with a multipolar +expansion, computes the in-cell capacitance C[0,j] (F/m) and inductance +L[0,j] (H/m) for a fixed reference conductor i=0, varying over the +conductors of interest listed in J_INDICES. + +Formulas (from InCellPotentials::getInCellCapacitanceUsingInnerRegion and +getInCellInductanceUsingInnerRegion in Results.cpp): + + C[0,j] = Q_j / (V_0|Vj=1 - _inner) * ε₀ + + L[0,j] = (A_0|Aj=1 - _inner) / I_j * μ₀ + +where: + Q_j = electric[j].ab[0][0] (monopole charge coefficient) + = electric[j].innerRegionAveragePotential + V_0 = electric[j].conductorPotentials[0] + + I_j = magnetic[j].ab[0][0] (monopole current coefficient) + = magnetic[j].innerRegionAveragePotential + A_0 = magnetic[j].conductorPotentials[0] + +Usage: + python compute_incell_parameters.py [CASE_DIR] + + CASE_DIR defaults to the realistic_case_fdtd_cell_centered_in_0 directory. +""" + +import json +import math +import os +import sys + +# Physical constants (match tulip/src/driver/constants.h) +EPSILON0_SI = 8.8541878176e-12 # F/m +MU0_SI = 4.0e-7 * math.pi # H/m + +# --------------------------------------------------------------------------- +# Configuration +# --------------------------------------------------------------------------- +DEFAULT_CASE_DIR = ( + "/home/luis/tulip/tmp_cases/realistic_case_fdtd_cell_centered_in_0" +) + +def load_json(case_dir: str) -> dict: + case_name = os.path.basename(case_dir.rstrip("/")) + json_path = os.path.join(case_dir, f"{case_name}.tulip.out.json") + with open(json_path) as fh: + return json.load(fh) + + +# --------------------------------------------------------------------------- +# Core computation (mirrors Results.cpp formulas) +# --------------------------------------------------------------------------- + +def capacitance(electric_solutions: list, i: int, j: int) -> float: + """C[i,j] in F/m.""" + sol_j = electric_solutions[j] + Q_j = sol_j["ab"][0][0] + avg_V_j = sol_j["innerRegionAveragePotential"] + V_i = sol_j["conductorPotentials"][i] + return Q_j / (V_i - avg_V_j) * EPSILON0_SI + + +def inductance(magnetic_solutions: list, i: int, j: int) -> float: + """L[i,j] in H/m.""" + sol_j = magnetic_solutions[j] + I_j = sol_j["ab"][0][0] + avg_A_j = sol_j["innerRegionAveragePotential"] + A_i = sol_j["conductorPotentials"][i] + return (A_i - avg_A_j) / I_j * MU0_SI + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + case_dir = sys.argv[1] if len(sys.argv) > 1 else DEFAULT_CASE_DIR + case_dir = os.path.abspath(case_dir) + + # Reference conductor + I_REF = 16 + + # Conductor indices j to evaluate + J_INDICES = [0, 16, 25, 30] + + + print(f"Case directory : {case_dir}") + print(f"Reference conductor i = {I_REF}") + print(f"Conductors j = {J_INDICES}") + + data = load_json(case_dir) + materials = data.get("materials", []) + mat_assoc = data.get("materialAssociations", []) + + if not materials: + print("No materials found in output JSON.") + return + + for material in materials: + mat_id = material["id"] + mat_type = material.get("type", "unknown") + mp = material.get("multipolarExpansion", {}) + + e_sols = mp.get("electric", []) + m_sols = mp.get("magnetic", []) + + if not e_sols or not m_sols: + print(f"[warn] Material {mat_id}: missing electric or magnetic solutions, skipping.") + continue + + assoc = next((a for a in mat_assoc if a["materialId"] == mat_id), None) + element_ids = assoc["elementIds"] if assoc else list(range(len(e_sols))) + + print(f"\n{'='*55}") + print(f" Material id={mat_id} type={mat_type}") + print(f" Total conductors: {len(element_ids)}") + 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}") + + for j in J_INDICES: + if j >= len(e_sols): + print(f" {j:>4} [index out of range]") + continue + C_val = capacitance(e_sols, I_REF, j) + L_val = inductance(m_sols, I_REF, j) + print(f" {j:>4} {C_val:>+18.6e} {L_val:>+18.6e}") + + print() + + +if __name__ == "__main__": + main() From 374a28727910204a907a29105bf773977d40da5d Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 13 May 2026 11:27:19 +0200 Subject: [PATCH 2/3] Makes sure that ids are used rather than positions --- scripts/compute_incell_parameters.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/scripts/compute_incell_parameters.py b/scripts/compute_incell_parameters.py index 1486ad1..400ed71 100644 --- a/scripts/compute_incell_parameters.py +++ b/scripts/compute_incell_parameters.py @@ -43,7 +43,7 @@ # Configuration # --------------------------------------------------------------------------- DEFAULT_CASE_DIR = ( - "/home/luis/tulip/tmp_cases/realistic_case_fdtd_cell_centered_in_0" + "/home/luis/workspace/tulip/tmp_cases/realistic_case_just_16_and_30" ) def load_json(case_dir: str) -> dict: @@ -75,6 +75,14 @@ def inductance(magnetic_solutions: list, i: int, j: int) -> float: 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: + return element_ids.index(element_id) + except ValueError: + return None + + # --------------------------------------------------------------------------- # Main # --------------------------------------------------------------------------- @@ -87,7 +95,7 @@ def main(): I_REF = 16 # Conductor indices j to evaluate - J_INDICES = [0, 16, 25, 30] + J_INDICES = [16, 30] print(f"Case directory : {case_dir}") @@ -117,6 +125,10 @@ def main(): assoc = next((a for a in mat_assoc if a["materialId"] == mat_id), None) element_ids = assoc["elementIds"] if assoc else list(range(len(e_sols))) + i_local = local_index_from_element_id(element_ids, I_REF) + if i_local is None: + print(f"[warn] Material {mat_id}: reference conductor id {I_REF} not in {element_ids}, skipping.") + continue print(f"\n{'='*55}") print(f" Material id={mat_id} type={mat_type}") print(f" Total conductors: {len(element_ids)}") @@ -125,11 +137,12 @@ def main(): print(f" {'-'*4} {'-'*18} {'-'*18}") for j in J_INDICES: - if j >= len(e_sols): - print(f" {j:>4} [index out of range]") + j_local = local_index_from_element_id(element_ids, j) + if j_local is None: + print(f" {j:>4} [id not in materialAssociations.elementIds]") continue - C_val = capacitance(e_sols, I_REF, j) - L_val = inductance(m_sols, I_REF, j) + C_val = capacitance(e_sols, i_local, j_local) + L_val = inductance(m_sols, i_local, j_local) print(f" {j:>4} {C_val:>+18.6e} {L_val:>+18.6e}") print() From bf43bc93bb50de55292f6c454872b7e27de943c6 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Thu, 14 May 2026 12:18:04 +0200 Subject: [PATCH 3/3] Addreses @adrianarce-elemwave comments --- scripts/compute_incell_parameters.py | 361 ++++++++++++++++++++------- 1 file changed, 269 insertions(+), 92 deletions(-) diff --git a/scripts/compute_incell_parameters.py b/scripts/compute_incell_parameters.py index 400ed71..13a0ac3 100644 --- a/scripts/compute_incell_parameters.py +++ b/scripts/compute_incell_parameters.py @@ -1,78 +1,71 @@ """ compute_incell_parameters.py -============================= +============================ Postprocessing script for a Tulip solved case. -Reads the .tulip.out.json file and, for each material with a multipolar -expansion, computes the in-cell capacitance C[0,j] (F/m) and inductance -L[0,j] (H/m) for a fixed reference conductor i=0, varying over the -conductors of interest listed in J_INDICES. +Reads a .tulip.out.json file and, for each material with a multipolar +expansion, computes the in-cell capacitance C[i,j] (F/m) and inductance +L[i,j] (H/m) for a fixed reference conductor i and all available conductors j +in that material association. Formulas (from InCellPotentials::getInCellCapacitanceUsingInnerRegion and getInCellInductanceUsingInnerRegion in Results.cpp): - C[0,j] = Q_j / (V_0|Vj=1 - _inner) * ε₀ + C[i,j] = Q_j / (V_i|Vj=1 - _inner) * epsilon0 - L[0,j] = (A_0|Aj=1 - _inner) / I_j * μ₀ + L[i,j] = (A_i|Aj=1 - _inner) / I_j * mu0 where: Q_j = electric[j].ab[0][0] (monopole charge coefficient) = electric[j].innerRegionAveragePotential - V_0 = electric[j].conductorPotentials[0] + V_i = electric[j].conductorPotentials[i] I_j = magnetic[j].ab[0][0] (monopole current coefficient) = magnetic[j].innerRegionAveragePotential - A_0 = magnetic[j].conductorPotentials[0] + A_i = magnetic[j].conductorPotentials[i] Usage: - python compute_incell_parameters.py [CASE_DIR] - - CASE_DIR defaults to the realistic_case_fdtd_cell_centered_in_0 directory. + python scripts/compute_incell_parameters.py + python scripts/compute_incell_parameters.py --run-tests """ +import argparse +import contextlib +import io import json import math import os import sys +import unittest # Physical constants (match tulip/src/driver/constants.h) -EPSILON0_SI = 8.8541878176e-12 # F/m -MU0_SI = 4.0e-7 * math.pi # H/m - -# --------------------------------------------------------------------------- -# Configuration -# --------------------------------------------------------------------------- -DEFAULT_CASE_DIR = ( - "/home/luis/workspace/tulip/tmp_cases/realistic_case_just_16_and_30" -) - -def load_json(case_dir: str) -> dict: - case_name = os.path.basename(case_dir.rstrip("/")) - json_path = os.path.join(case_dir, f"{case_name}.tulip.out.json") - with open(json_path) as fh: - return json.load(fh) +EPSILON0_SI = 8.8541878176e-12 # F/m +MU0_SI = 4.0e-7 * math.pi # H/m -# --------------------------------------------------------------------------- -# Core computation (mirrors Results.cpp formulas) -# --------------------------------------------------------------------------- +def load_json(json_path: str) -> dict: + """Load a Tulip output JSON file from disk.""" + abs_path = os.path.abspath(json_path) + with open(abs_path, encoding="utf-8") as fh: + return json.load(fh) + def capacitance(electric_solutions: list, i: int, j: int) -> float: - """C[i,j] in F/m.""" - sol_j = electric_solutions[j] - Q_j = sol_j["ab"][0][0] - avg_V_j = sol_j["innerRegionAveragePotential"] - V_i = sol_j["conductorPotentials"][i] - return Q_j / (V_i - avg_V_j) * EPSILON0_SI + """Compute C[i,j] in F/m.""" + sol_j = electric_solutions[j] + q_j = sol_j["ab"][0][0] + avg_v_j = sol_j["innerRegionAveragePotential"] + v_i = sol_j["conductorPotentials"][i] + return q_j / (v_i - avg_v_j) * EPSILON0_SI def inductance(magnetic_solutions: list, i: int, j: int) -> float: - """L[i,j] in H/m.""" - sol_j = magnetic_solutions[j] - I_j = sol_j["ab"][0][0] - avg_A_j = sol_j["innerRegionAveragePotential"] - A_i = sol_j["conductorPotentials"][i] - return (A_i - avg_A_j) / I_j * MU0_SI + """Compute L[i,j] in H/m.""" + sol_j = magnetic_solutions[j] + i_j = sol_j["ab"][0][0] + avg_a_j = sol_j["innerRegionAveragePotential"] + 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): @@ -83,70 +76,254 @@ def local_index_from_element_id(element_ids: list, element_id: int): return None -# --------------------------------------------------------------------------- -# Main -# --------------------------------------------------------------------------- +def element_ids_for_material(material_id: int, mat_assoc: list, fallback_size: int) -> list: + """Return material-associated conductor ids, or fallback to local indices.""" + assoc = next((a for a in mat_assoc if a.get("materialId") == material_id), None) + if assoc is not None and "elementIds" in assoc: + return assoc["elementIds"] + return list(range(fallback_size)) -def main(): - case_dir = sys.argv[1] if len(sys.argv) > 1 else DEFAULT_CASE_DIR - case_dir = os.path.abspath(case_dir) - # Reference conductor - I_REF = 16 +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.""" + 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", []) - # Conductor indices j to evaluate - J_INDICES = [16, 30] + if not e_sols or not m_sols: + return None, "missing electric or 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: + return None, f"reference conductor id {i_ref} not in {element_ids}" - print(f"Case directory : {case_dir}") - print(f"Reference conductor i = {I_REF}") - print(f"Conductors j = {J_INDICES}") - - data = load_json(case_dir) + rows = [] + for j in element_ids: + 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) + rows.append((j, c_val, l_val)) + + return { + "mat_id": mat_id, + "mat_type": mat_type, + "element_ids": element_ids, + "rows": rows, + }, None + + +def print_material_report(report: dict, i_ref: int): + """Print computed rows for one material.""" + 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"{'=' * 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}") + + for j, c_val, l_val in report["rows"]: + print(f" {j:>4} {c_val:>+18.6e} {l_val:>+18.6e}") + + +def parse_args(argv: list): + """Parse command-line arguments.""" + parser = argparse.ArgumentParser( + description="Compute in-cell C/L parameters from a .tulip.out.json file." + ) + parser.add_argument( + "--run-tests", + action="store_true", + help="Run the unit tests embedded in this file.", + ) + parser.add_argument( + "json_path", + nargs="?", + help="Path to a .tulip.out.json file.", + ) + parser.add_argument( + "i_ref", + nargs="?", + type=int, + help="Reference conductor id i.", + ) + + args = parser.parse_args(argv) + if args.run_tests: + return args + + 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") + return args + + +def run_self_tests() -> int: + """Execute same-file unit tests.""" + test_program = unittest.main(module=__name__, argv=[sys.argv[0]], exit=False) + return 0 if test_program.result.wasSuccessful() else 1 + + +def run(argv=None) -> int: + """Entrypoint used by both CLI and tests.""" + args = parse_args(sys.argv[1:] if argv is None else argv) + if args.run_tests: + return run_self_tests() + + json_path = os.path.abspath(args.json_path) + i_ref = args.i_ref + + print(f"Output JSON path : {json_path}") + print(f"Reference conductor i : {i_ref}") + + data = load_json(json_path) materials = data.get("materials", []) mat_assoc = data.get("materialAssociations", []) if not materials: print("No materials found in output JSON.") - return + return 0 for material in materials: - mat_id = material["id"] - mat_type = material.get("type", "unknown") - mp = material.get("multipolarExpansion", {}) - - e_sols = mp.get("electric", []) - m_sols = mp.get("magnetic", []) - - if not e_sols or not m_sols: - print(f"[warn] Material {mat_id}: missing electric or magnetic solutions, skipping.") - continue - - assoc = next((a for a in mat_assoc if a["materialId"] == mat_id), None) - element_ids = assoc["elementIds"] if assoc else list(range(len(e_sols))) - - i_local = local_index_from_element_id(element_ids, I_REF) - if i_local is None: - print(f"[warn] Material {mat_id}: reference conductor id {I_REF} not in {element_ids}, skipping.") + report, warn = compute_material_rows(material, mat_assoc, i_ref) + if warn is not None: + print(f"[warn] Material {material.get('id')}: {warn}, skipping.") continue - print(f"\n{'='*55}") - print(f" Material id={mat_id} type={mat_type}") - print(f" Total conductors: {len(element_ids)}") - 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}") - - for j in J_INDICES: - j_local = local_index_from_element_id(element_ids, j) - if j_local is None: - print(f" {j:>4} [id not in materialAssociations.elementIds]") - continue - C_val = capacitance(e_sols, i_local, j_local) - L_val = inductance(m_sols, i_local, j_local) - print(f" {j:>4} {C_val:>+18.6e} {L_val:>+18.6e}") - + print_material_report(report, i_ref) print() + return 0 + + +class TestComputeInCellParameters(unittest.TestCase): + """Unit tests for compute_incell_parameters.py.""" + + def test_local_index_lookup(self): + element_ids = [16, 30, 42] + self.assertEqual(local_index_from_element_id(element_ids, 30), 1) + self.assertIsNone(local_index_from_element_id(element_ids, 99)) + + def test_compute_material_rows_uses_all_conductors(self): + material = { + "id": 7, + "type": "dielectric", + "multipolarExpansion": { + "electric": [ + { + "ab": [[1.0]], + "innerRegionAveragePotential": 0.0, + "conductorPotentials": [2.0, 3.0], + }, + { + "ab": [[2.0]], + "innerRegionAveragePotential": 1.0, + "conductorPotentials": [3.0, 4.0], + }, + ], + "magnetic": [ + { + "ab": [[1.0]], + "innerRegionAveragePotential": 0.0, + "conductorPotentials": [2.0, 3.0], + }, + { + "ab": [[2.0]], + "innerRegionAveragePotential": 1.0, + "conductorPotentials": [3.0, 4.0], + }, + ], + }, + } + mat_assoc = [{"materialId": 7, "elementIds": [16, 30]}] + + report, warn = compute_material_rows(material, mat_assoc, 16) + self.assertIsNone(warn) + self.assertIsNotNone(report) + self.assertEqual([row[0] for row in report["rows"]], [16, 30]) + + def test_reference_case_values_from_embedded_fixture(self): + material = { + "id": 1, + "type": "unshieldedMultiwire", + "multipolarExpansion": { + "electric": [ + { + "ab": [ + [0.8419073593680602, 0.0], + [4.7094817324223296e-07, -7.437578541169077e-07], + [-2.254972536944232e-06, -1.3845721768518446e-06], + [-1.7502712357816652e-09, 1.6065198400637674e-09], + ], + "conductorPotentials": [1.0, 0.7525321728251249], + "expansionCenter": [0.0017138274803646352, -0.0027066091802577024], + "innerRegionAveragePotential": 0.6100537777863153, + }, + { + "ab": [ + [1.0308546965043754, 0.0], + [4.236110742075576e-07, -8.85263389452124e-07], + [-1.301845581534472e-07, -8.009832224471683e-08], + [3.712118282384746e-10, -3.4152957909780026e-10], + ], + "conductorPotentials": [0.9213984270827741, 1.0], + "expansionCenter": [0.0023197106176927533, -0.004847736541849076], + "innerRegionAveragePotential": 0.7447102973543827, + }, + ], + "magnetic": [ + { + "ab": [ + [0.8419073593680602, 0.0], + [4.7094817324223296e-07, -7.437578541169077e-07], + [-2.254972536944232e-06, -1.3845721768518446e-06], + [-1.7502712357816652e-09, 1.6065198400637674e-09], + ], + "conductorPotentials": [1.0, 0.7525321728251249], + "expansionCenter": [0.0017138274803646352, -0.0027066091802577024], + "innerRegionAveragePotential": 0.6100537777863153, + }, + { + "ab": [ + [1.0308546965043754, 0.0], + [4.236110742075576e-07, -8.85263389452124e-07], + [-1.301845581534472e-07, -8.009832224471683e-08], + [3.712118282384746e-10, -3.4152957909780026e-10], + ], + "conductorPotentials": [0.9213984270827741, 1.0], + "expansionCenter": [0.0023197106176927533, -0.004847736541849076], + "innerRegionAveragePotential": 0.7447102973543827, + }, + ], + }, + } + mat_assoc = [{"elementIds": [16, 30], "materialId": 1}] + + report, warn = compute_material_rows(material, mat_assoc, 16) + self.assertIsNone(warn) + self.assertIsNotNone(report) + + rows = {j: (c_val, l_val) for j, c_val, l_val in report["rows"]} + self.assertIn(16, rows) + self.assertIn(30, rows) + + self.assertAlmostEqual(rows[16][0], 1.9116497250688996e-11, delta=1e-20) + self.assertAlmostEqual(rows[16][1], 5.820365736777192e-07, delta=1e-15) + 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): + with contextlib.redirect_stderr(io.StringIO()): + with self.assertRaises(SystemExit): + parse_args([]) + + def test_parse_args_happy_path(self): + args = parse_args(["case.tulip.out.json", "16"]) + self.assertEqual(args.json_path, "case.tulip.out.json") + self.assertEqual(args.i_ref, 16) + if __name__ == "__main__": - main() + sys.exit(run())