From 02f43d1545b72be2f52c9b16f18f17471841506d Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Wed, 18 Feb 2026 16:16:07 +0100 Subject: [PATCH 01/10] WIP --- src_json_parser/parser_tools.F90 | 7 +++ src_json_parser/smbjson.F90 | 81 ++++++++++++++++++++++-------- src_mtln/mtl.F90 | 42 ++++++++++++++-- src_mtln/mtln_types.F90 | 2 + src_mtln/preprocess.F90 | 2 +- test/mtln/mtln_testingTools.F90 | 2 +- test/mtln/test_mtl.F90 | 2 +- test/mtln/test_mtl_bundle.F90 | 2 +- test/pyWrapper/test_full_system.py | 46 +++++++++-------- test/pyWrapper/utils.py | 2 +- 10 files changed, 137 insertions(+), 51 deletions(-) diff --git a/src_json_parser/parser_tools.F90 b/src_json_parser/parser_tools.F90 index b8dc2a1d..3c6174c0 100644 --- a/src_json_parser/parser_tools.F90 +++ b/src_json_parser/parser_tools.F90 @@ -191,6 +191,13 @@ function getPixelFromElementId(mesh, id) result(res) end function + function vectorSize(vector) result(res) + real, dimension(:), intent(in) :: vector + integer :: n, res + n = size(vector,1) + res = n + end function + function vectorToDiagonalMatrix(vector) result(res) real, dimension(:), intent(in) :: vector real, dimension(:, :), allocatable :: res diff --git a/src_json_parser/smbjson.F90 b/src_json_parser/smbjson.F90 index b13d5e11..3f3a3357 100644 --- a/src_json_parser/smbjson.F90 +++ b/src_json_parser/smbjson.F90 @@ -1225,8 +1225,13 @@ function readMoreProbes(this) result (res) type(json_value_ptr), dimension(:), allocatable :: ps integer :: i +#ifdef CompileWithMTLN + character (len=*), dimension(2), parameter :: validTypes = & + [J_PR_TYPE_POINT, J_PR_TYPE_LINE] +#else character (len=*), dimension(3), parameter :: validTypes = & [J_PR_TYPE_POINT, J_PR_TYPE_WIRE, J_PR_TYPE_LINE] +#endif logical :: found character (len=:), allocatable :: fieldLbl, probeLbl integer :: filtered_size, n @@ -2518,18 +2523,14 @@ function readMTLN(this) result (mtln_res) class(parser_t) :: this type(mtln_t) :: mtln_res type(fhash_tbl_t) :: elemIdToPosition, elemIdToCable, connIdToConnector - type(materialAssociation_t), dimension(:), allocatable :: wires, cables + type(materialAssociation_t), dimension(:), allocatable :: cables class(cable_t), pointer :: ptr, read_cable integer :: i - wires = this%getMaterialAssociations([J_MAT_TYPE_WIRE]) - if (size(wires) /=0 ) then - call WarnErrReport("ERROR: material type 'wire' is not allowed if compiled with MTLN", .true.) - end if - cables = this%getMaterialAssociations([ & J_MAT_TYPE_SHIELDED_MULTIWIRE//' ',& - J_MAT_TYPE_UNSHIELDED_MULTIWIRE ]) + J_MAT_TYPE_UNSHIELDED_MULTIWIRE ,& + J_MAT_TYPE_WIRE//' ' ]) ! spaces are needed to make strings have same length. ! Why? Because of FORTRAN! It only accepts fixed length strings for arrays. @@ -2589,6 +2590,8 @@ function assignParentCable(cable, cables) result(res) end if case (J_MAT_TYPE_UNSHIELDED_MULTIWIRE) res => null() + case (J_MAT_TYPE_WIRE) + res => null() case default call WarnErrReport('ERROR: Material type not recognized', .true.) end select @@ -2612,6 +2615,8 @@ function assignConductorInParent(cable) result(res) end if case (J_MAT_TYPE_UNSHIELDED_MULTIWIRE) res = 0 + case (J_MAT_TYPE_WIRE) + res = 0 case default call WarnErrReport('ERROR: Material type not recognized', .true.) end select @@ -2712,7 +2717,8 @@ function buildNetworks() result(res) allocate(networks_coordinates(0)) ! cables = [ this%getMaterialAssociations([J_MAT_TYPE_WIRE]), & cables = [this%getMaterialAssociations([J_MAT_TYPE_UNSHIELDED_MULTIWIRE]), & - this%getMaterialAssociations([J_MAT_TYPE_SHIELDED_MULTIWIRE]) ] + this%getMaterialAssociations([J_MAT_TYPE_SHIELDED_MULTIWIRE]) ,& + this%getMaterialAssociations([J_MAT_TYPE_WIRE]) ] do i = 1, size(cables) elemIds = cables(i)%elementIds terminations_ini => getTerminationsOnSide(cables(i)%initialTerminalId) @@ -3264,7 +3270,8 @@ function countOutputProbes(probes) result(res) mAs = this%getMaterialAssociations([ & J_MAT_TYPE_SHIELDED_MULTIWIRE//' ',& - J_MAT_TYPE_UNSHIELDED_MULTIWIRE ]) + J_MAT_TYPE_UNSHIELDED_MULTIWIRE ,& + J_MAT_TYPE_WIRE//' ' ]) do i = 1, size(mAs) polyline = this%mesh%getPolyline(mAs(i)%elementIds(1)) @@ -3303,7 +3310,8 @@ logical function isProbeDefinedOnMultiwire(p) mAs = this%getMaterialAssociations([ & J_MAT_TYPE_SHIELDED_MULTIWIRE//' ',& - J_MAT_TYPE_UNSHIELDED_MULTIWIRE ]) + J_MAT_TYPE_UNSHIELDED_MULTIWIRE ,& + J_MAT_TYPE_WIRE//' ' ]) do i = 1, size(mAs) polyline = this%mesh%getPolyline(mAs(i)%elementIds(1)) @@ -3336,7 +3344,8 @@ function getPolylineElemIdOfMultiwireProbe(p) result(res) mAs = this%getMaterialAssociations([ & J_MAT_TYPE_SHIELDED_MULTIWIRE//' ',& - J_MAT_TYPE_UNSHIELDED_MULTIWIRE ]) + J_MAT_TYPE_UNSHIELDED_MULTIWIRE ,& + J_MAT_TYPE_WIRE//' ' ]) allocate(res(0)) do i = 1, size(mAs) polyline = this%mesh%getPolyline(mAs(i)%elementIds(1)) @@ -3478,7 +3487,7 @@ function readMTLNCable(j_cable, despl) result(res) res%transfer_impedance = buildTransferImpedance(material) call assignPULProperties(res, material, size(j_cable%elementIds)) end select - case (J_MAT_TYPE_UNSHIELDED_MULTIWIRE) + case (J_MAT_TYPE_UNSHIELDED_MULTIWIRE, J_MAT_TYPE_WIRE) allocate(unshielded_multiwire_t :: res) select type(res) type is(unshielded_multiwire_t) @@ -3548,10 +3557,12 @@ subroutine assignInCellProperties(res, mat, n) type(json_value_ptr) :: mat type(json_value), pointer :: multipolarExpansionPtr integer, intent(in) :: n + integer :: m real, dimension(:,:), allocatable :: null_matrix logical :: found logical :: areFixedInCell logical :: areMultipolarInCell + logical :: hasRadius allocate(null_matrix(n,n), source = 0.0) @@ -3560,34 +3571,53 @@ subroutine assignInCellProperties(res, mat, n) this%existsAt(mat%p, J_MAT_MULTIWIRE_CAPACITANCE) areMultipolarInCell = & this%existsAt(mat%p, J_MAT_MULTIWIRE_MULTIPOLAR_EXPANSION) - - if ((areFixedInCell .and. areMultipolarInCell) .or. & - (.not. areFixedInCell .and. .not. areMultipolarInCell) ) then - call WarnErrReport( & - "Unshielded multiwires in cell properties must be defined by fixed OR multipolarExpansions, but not both.", .true.) + hasRadius = & + this%existsAt(mat%p, J_MAT_WIRE_RADIUS) .and. & + this%getRealAt(mat%p, J_MAT_WIRE_RADIUS, default = 0.0) /= 0 + + if (.not. hasRadius) then + if ((areFixedInCell .and. areMultipolarInCell) .or. & + (.not. areFixedInCell .and. .not. areMultipolarInCell) ) then + call WarnErrReport( & + "Unshielded multiwires in cell properties must be defined by fixed OR multipolarExpansions, but not both.", .true.) + end if end if - if (areFixedInCell) then res%cell_inductance_per_meter = this%getMatrixAt(mat%p, J_MAT_MULTIWIRE_INDUCTANCE,found) res%cell_capacitance_per_meter = this%getMatrixAt(mat%p, J_MAT_MULTIWIRE_CAPACITANCE,found) allocate(res%multipolar_expansion(0)) - else + else if (areMultipolarInCell) then res%cell_inductance_per_meter = null_matrix res%cell_capacitance_per_meter = null_matrix call this%core%get(mat%p, J_MAT_MULTIWIRE_MULTIPOLAR_EXPANSION, multipolarExpansionPtr) allocate(res%multipolar_expansion(1)) res%multipolar_expansion(1) = readMultipolarExpansion(multipolarExpansionPtr) + else if (hasRadius) then + res%cell_inductance_per_meter = null_matrix + res%cell_capacitance_per_meter = null_matrix + res%radius = this%getRealAt(mat%p, J_MAT_WIRE_RADIUS, default = 0.0) end if + if (this%existsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE)) then - res%resistance_per_meter = & - vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found)) + m = vectorSize(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found)) + ! resistance = ... + if (m == n) then + res%resistance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found)) + else + res%resistance_per_meter = null_matrix + end if else res%resistance_per_meter = null_matrix end if if (this%existsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE)) then - res%conductance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found)) + m = vectorSize(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found)) + if (m == n) then + res%conductance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found)) + else + res%conductance_per_meter = null_matrix + end if else res%conductance_per_meter = null_matrix end if @@ -3683,10 +3713,17 @@ function buildSegments(j_cable, despl) result(res) select case(abs(res(i)%orientation)) case(DIR_X) res(i)%dualBox = getdualBoxYZ(res(i), despl) + res(i)%d1 = despl%desY(res(i)%y) + res(i)%d2 = despl%desZ(res(i)%z) case(DIR_Y) res(i)%dualBox = getdualBoxZX(res(i), despl) + res(i)%d1 = despl%desZ(res(i)%z) + res(i)%d1 = despl%desX(res(i)%x) case(DIR_Z) res(i)%dualBox = getdualBoxXY(res(i), despl) + res(i)%d1 = despl%desX(res(i)%x) + res(i)%d2 = despl%desY(res(i)%y) + end select end if prevOr = abs(res(i)%orientation) diff --git a/src_mtln/mtl.F90 b/src_mtln/mtl.F90 index 4083848f..df6ee8d1 100644 --- a/src_mtln/mtl.F90 +++ b/src_mtln/mtl.F90 @@ -6,7 +6,9 @@ module mtl_mod use mtln_types_mod, only: segment_t, multipolar_expansion_t use multipolar_expansion_mod, only: getCellCapacitanceOnBox, getCellInductanceOnBox #ifdef CompileWithMPI - use fdetypes, only: SUBCOMM_MPI, REALSIZE, INTEGERSIZE + use fdetypes, only: SUBCOMM_MPI, REALSIZE, INTEGERSIZE, pi, mu_vacuum, RKIND_wires +#else + use fdetypes, only: pi, mu_vacuum, RKIND_wires #endif implicit none #ifdef CompileWithMPI @@ -59,6 +61,7 @@ module mtl_mod procedure :: checkTimeStep procedure :: allocatePULMatrices procedure :: computeLCParameters + procedure :: computeLCParametersFromRadius procedure :: initLC procedure :: initRG procedure :: initDirections @@ -163,7 +166,7 @@ function mtl_shielded(lpul, cpul, rpul, gpul, & function mtl_unshielded(lpul, cpul, rpul, gpul, & step_size, name, segments, & - dt, multipolar_expansion, & + dt, multipolar_expansion, radius, & layer_indices, bundle_in_layer, alloc_z) result(res) real, intent(in), dimension(:,:) :: lpul, cpul, rpul, gpul real, intent(in), dimension(:) :: step_size @@ -171,7 +174,7 @@ function mtl_unshielded(lpul, cpul, rpul, gpul, & type(segment_t), dimension(:), allocatable, intent(in) :: segments real, intent(in) :: dt type(multipolar_expansion_t), dimension(:), allocatable :: multipolar_expansion - + real, intent(in) :: radius integer (kind=4), allocatable, dimension(:,:), intent(in), optional :: layer_indices logical, optional :: bundle_in_layer integer(kind=4), dimension (2), intent(in), optional :: alloc_z @@ -203,7 +206,9 @@ function mtl_unshielded(lpul, cpul, rpul, gpul, & call res%allocatePULMatrices() if (size(multipolar_expansion) /= 0) then call res%computeLCParameters(multipolar_expansion(1)) - else + else if (radius /= 0.0) then + call res%computeLCParametersFromRadius(radius) + else call res%initLC(lpul, cpul) end if call res%initRG(rpul, gpul) @@ -261,6 +266,35 @@ subroutine computeLCParameters(this, multipolar_expansion) this%cpul(size(this%segments)+1, :,:) = this%cpul(size(this%segments), :,:) end subroutine + subroutine computeLCParametersFromRadius(this, rad) + class(mtl_t) :: this + real, intent(in) :: rad + REAL (KIND=RKIND_wires) :: invMu + integer :: i + real :: d1, d2 + invMu = 1.0/mu_vacuum + do i = 1, size(this%segments) + d1 = this%segments(i)%d1 + d2 = this%segments(i)%d2 + this%lpul(i,:,:) = & + (1.0_RKIND_wires / (4.0_RKIND_wires * pi*invMu))*(log((d1**2.0_RKIND_wires +d2**2.0_RKIND_wires )/(4.0_RKIND_wires *rad**2.0_RKIND_wires ))+ & + d1/d2*atan(d2/d1) + & + d2/d1*atan(d1/d2) + pi*rad**2.0_RKIND_wires /(d2*d1)-3.0_RKIND_wires) + + if ((rad < 0.3_RKIND_wires *d1).or.(rad < 0.3_RKIND_wires *d2)) then + this%lpul(i,:,:) = this%lpul(i,:,:) - 0.57_RKIND_wires/(4.0_RKIND_wires * pi*invMu) + endif + + if ((rad > 0.3_RKIND_wires *d1).or.(rad > 0.3_RKIND_wires *d2)) then + this%lpul(i,:,:) = this%lpul(i,:,:) & + /(1.0_RKIND_wires-pi*rad**2.0_RKIND_wires /(d1*d2)) + endif + this%cpul(i,:,:) = 1.0/this%lpul(i,:,:) + enddo + this%cpul(size(this%segments)+1, :,:) = this%cpul(size(this%segments), :,:) + + end subroutine + subroutine initDirections(this) class(mtl_t) :: this diff --git a/src_mtln/mtln_types.F90 b/src_mtln/mtln_types.F90 index 3d4ade17..4672cdd1 100644 --- a/src_mtln/mtln_types.F90 +++ b/src_mtln/mtln_types.F90 @@ -180,6 +180,7 @@ module mtln_types_mod type, extends(direction_t) :: segment_t type(box_2d_t) :: dualBox + real :: d1, d2 end type @@ -200,6 +201,7 @@ module mtln_types_mod real, allocatable, dimension(:,:) :: cell_capacitance_per_meter real, allocatable, dimension(:,:) :: resistance_per_meter real, allocatable, dimension(:,:) :: conductance_per_meter + real :: radius = 0.0 ! should multipolar expansion be always present, instead of allocatable, ! but check if it is being used using the size of the field_reconstruction? type(multipolar_expansion_t), dimension(:), allocatable :: multipolar_expansion diff --git a/src_mtln/preprocess.F90 b/src_mtln/preprocess.F90 index 93fc0920..8453fe1c 100644 --- a/src_mtln/preprocess.F90 +++ b/src_mtln/preprocess.F90 @@ -274,7 +274,7 @@ function buildLineFromCable(cable, dt, layer_indices, bundle_in_layer, alloc_z) res = mtl_unshielded( lpul = cable%cell_inductance_per_meter, cpul = cable%cell_capacitance_per_meter, & rpul = cable%resistance_per_meter, gpul = cable%conductance_per_meter, & step_size = cable%step_size, name = cable%name, segments = cable%segments,& - dt = dt, multipolar_expansion = cable%multipolar_expansion & + dt = dt, multipolar_expansion = cable%multipolar_expansion, radius = cable%radius & #ifdef CompileWithMPI ,layer_indices = layer_indices, bundle_in_layer = bundle_in_layer, alloc_z = alloc_z & #endif diff --git a/test/mtln/mtln_testingTools.F90 b/test/mtln/mtln_testingTools.F90 index 93233afa..1567f48d 100644 --- a/test/mtln/mtln_testingTools.F90 +++ b/test/mtln/mtln_testingTools.F90 @@ -82,7 +82,7 @@ type(mtl_t) function buildLineWithNConductors(n,name, parent_name, conductor_in_ end if res = mtl_shielded(lpul, cpul, rpul, gpul, step_size, name, segments, time_step, parent, conductor, Zt) else if (type == MTL_TYPE_UNSHIELDED) then - res = mtl_unshielded(lpul, cpul, rpul, gpul, step_size, name, segments, time_step, mE) + res = mtl_unshielded(lpul, cpul, rpul, gpul, step_size, name, segments, time_step, mE, radius = 0.0) else write(*,*) 'Unrecognized line type' end if diff --git a/test/mtln/test_mtl.F90 b/test/mtln/test_mtl.F90 index 5bbe470c..5b6fab90 100644 --- a/test/mtln/test_mtl.F90 +++ b/test/mtln/test_mtl.F90 @@ -58,7 +58,7 @@ integer function test_mtl_init_homogeneous() bind(C) result(error_cnt) call comparePULMatrices(error_cnt, line%cpul, cpul) call comparePULMatrices(error_cnt, line%rpul, rpul) call comparePULMatrices(error_cnt, line%gpul, gpul) - line = mtl_unshielded(lpul, cpul, rpul, gpul, step_size, name, segments=segments, dt = 1e-12, multipolar_expansion = mE) + line = mtl_unshielded(lpul, cpul, rpul, gpul, step_size, name, segments=segments, dt = 1e-12, multipolar_expansion = mE, radius = 0.0) call comparePULMatrices(error_cnt, line%lpul, lpul) call comparePULMatrices(error_cnt, line%cpul, cpul) call comparePULMatrices(error_cnt, line%rpul, rpul) diff --git a/test/mtln/test_mtl_bundle.F90 b/test/mtln/test_mtl_bundle.F90 index 1d9c10bc..d6ffd6b4 100644 --- a/test/mtln/test_mtl_bundle.F90 +++ b/test/mtln/test_mtl_bundle.F90 @@ -36,7 +36,7 @@ integer function test_mtl_bundle_init() bind(C) result(error_cnt) parent_name = "line_out", conductor_in_parent = 1, & transfer_impedance = Zt) mtl_out = mtl_unshielded(l1, c1, r1, g1, step_size, name = "line_out", segments = segments, dt = 1e-11, & - multipolar_expansion = mE ) + multipolar_expansion = mE, radius = 0.0 ) allocate(levels(1)%lines(1)) allocate(levels(2)%lines(1)) diff --git a/test/pyWrapper/test_full_system.py b/test/pyWrapper/test_full_system.py index e3ca8827..5ecae5a3 100644 --- a/test/pyWrapper/test_full_system.py +++ b/test/pyWrapper/test_full_system.py @@ -174,22 +174,27 @@ def test_coated_antenna(tmp_path): p_solved['current_0'].to_numpy(), rtol=0.0, atol=10e-8) + def test_holland(tmp_path): fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][0] = createWire(id = 1, r = 0.02) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) + # if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): + # solver['materials'][0] = createWire(id = 1, r = 0.02) + # elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): + # solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) + solver['materials'][0] = createWire(id = 1, r = 0.02) solver.run() + p_wire = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) + + solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) + solver.cleanUp() + solver.run() + p_unshielded = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) - probe_current = solver.getSolvedProbeFilenames("mid_point")[0] - probe_files = [probe_current] - p_solved = Probe(probe_files[0]) expected_f = json.load(open(OUTPUTS_FOLDER+'holland1981.fdtd_mid_point_Wz_11_11_12_s2_12.json')) expected_t, expected_i = np.array([]), np.array([]) @@ -197,19 +202,20 @@ def test_holland(tmp_path): expected_t = np.append(expected_t, float(data['value'][0])) expected_i = np.append(expected_i, float(data['value'][1])) - expected_i_interp = np.interp(p_solved['time']-3.05*1e-9, expected_t, expected_i) - - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - assert np.allclose( - expected_i_interp, - -p_solved['current'], - rtol=1e-4, atol=5e-5) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - assert np.allclose( - expected_i_interp, - p_solved['current_0'], - rtol=1e-4, atol=5e-5) - + expected_i_interp = np.interp(p_wire['time']-3.05*1e-9, expected_t, expected_i) + + # if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): + # assert np.allclose( + # expected_i_interp, + # -p_solved['current'], + # rtol=1e-4, atol=5e-5) + # elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): + # assert np.allclose( + # expected_i_interp, + # p_solved['current_0'], + # rtol=1e-4, atol=5e-5) + assert np.allclose(expected_i_interp, p_wire['current_0'], rtol=1e-4, atol=5e-5) + assert np.allclose(expected_i_interp, p_unshielded['current_0'], rtol=1e-4, atol=5e-5) @no_mtln_skip diff --git a/test/pyWrapper/utils.py b/test/pyWrapper/utils.py index 7f1d1857..b86eb333 100644 --- a/test/pyWrapper/utils.py +++ b/test/pyWrapper/utils.py @@ -34,7 +34,7 @@ # Use of absolute path to avoid conflicts when changing directory. if platform == "linux": - SEMBA_EXE = os.path.join(os.getcwd(), 'build', 'bin', 'semba-fdtd') + SEMBA_EXE = os.path.join(os.getcwd(), 'build-dbg', 'bin', 'semba-fdtd') elif platform == "win32": SEMBA_EXE = os.path.join(os.getcwd(), 'build', 'bin', 'semba-fdtd.exe') From 00c219b13003067a26b16c6add639d34d0b6fa00 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 09:39:44 +0100 Subject: [PATCH 02/10] Bouyated inductance model implemented for wires in mtln --- src_json_parser/parser_tools.F90 | 6 ------ src_json_parser/smbjson.F90 | 15 +++++++++------ src_mtln/mtl.F90 | 9 +++++---- 3 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src_json_parser/parser_tools.F90 b/src_json_parser/parser_tools.F90 index 3c6174c0..e2f6a737 100644 --- a/src_json_parser/parser_tools.F90 +++ b/src_json_parser/parser_tools.F90 @@ -191,12 +191,6 @@ function getPixelFromElementId(mesh, id) result(res) end function - function vectorSize(vector) result(res) - real, dimension(:), intent(in) :: vector - integer :: n, res - n = size(vector,1) - res = n - end function function vectorToDiagonalMatrix(vector) result(res) real, dimension(:), intent(in) :: vector diff --git a/src_json_parser/smbjson.F90 b/src_json_parser/smbjson.F90 index 3f3a3357..b067cac1 100644 --- a/src_json_parser/smbjson.F90 +++ b/src_json_parser/smbjson.F90 @@ -3563,7 +3563,8 @@ subroutine assignInCellProperties(res, mat, n) logical :: areFixedInCell logical :: areMultipolarInCell logical :: hasRadius - + real, dimension(:), allocatable :: r, c + allocate(null_matrix(n,n), source = 0.0) areFixedInCell = & @@ -3596,13 +3597,13 @@ subroutine assignInCellProperties(res, mat, n) else if (hasRadius) then res%cell_inductance_per_meter = null_matrix res%cell_capacitance_per_meter = null_matrix + allocate(res%multipolar_expansion(0)) res%radius = this%getRealAt(mat%p, J_MAT_WIRE_RADIUS, default = 0.0) end if if (this%existsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE)) then - m = vectorSize(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found)) - ! resistance = ... - if (m == n) then + r = this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found) + if (found) then res%resistance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found)) else res%resistance_per_meter = null_matrix @@ -3612,8 +3613,8 @@ subroutine assignInCellProperties(res, mat, n) end if if (this%existsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE)) then - m = vectorSize(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found)) - if (m == n) then + c = this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found) + if (found) then res%conductance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found)) else res%conductance_per_meter = null_matrix @@ -3709,6 +3710,8 @@ function buildSegments(j_cable, despl) result(res) res(i)%orientation = linels(i)%orientation if (prevOr == abs(res(i)%orientation)) then res(i)%dualBox = res(i-1)%dualBox + res(i)%d1 = res(i-1)%d1 + res(i)%d2 = res(i-1)%d2 else select case(abs(res(i)%orientation)) case(DIR_X) diff --git a/src_mtln/mtl.F90 b/src_mtln/mtl.F90 index df6ee8d1..b4883e08 100644 --- a/src_mtln/mtl.F90 +++ b/src_mtln/mtl.F90 @@ -6,9 +6,9 @@ module mtl_mod use mtln_types_mod, only: segment_t, multipolar_expansion_t use multipolar_expansion_mod, only: getCellCapacitanceOnBox, getCellInductanceOnBox #ifdef CompileWithMPI - use fdetypes, only: SUBCOMM_MPI, REALSIZE, INTEGERSIZE, pi, mu_vacuum, RKIND_wires + use fdetypes, only: SUBCOMM_MPI, REALSIZE, INTEGERSIZE, pi, mu_vacuum, epsilon_vacuum, RKIND_wires #else - use fdetypes, only: pi, mu_vacuum, RKIND_wires + use fdetypes, only: pi, mu_vacuum, epsilon_vacuum, RKIND_wires #endif implicit none #ifdef CompileWithMPI @@ -269,10 +269,11 @@ subroutine computeLCParameters(this, multipolar_expansion) subroutine computeLCParametersFromRadius(this, rad) class(mtl_t) :: this real, intent(in) :: rad - REAL (KIND=RKIND_wires) :: invMu + REAL (KIND=RKIND_wires) :: invMu, c integer :: i real :: d1, d2 invMu = 1.0/mu_vacuum + c = 1/sqrt(epsilon_vacuum*mu_vacuum) do i = 1, size(this%segments) d1 = this%segments(i)%d1 d2 = this%segments(i)%d2 @@ -289,7 +290,7 @@ subroutine computeLCParametersFromRadius(this, rad) this%lpul(i,:,:) = this%lpul(i,:,:) & /(1.0_RKIND_wires-pi*rad**2.0_RKIND_wires /(d1*d2)) endif - this%cpul(i,:,:) = 1.0/this%lpul(i,:,:) + this%cpul(i,:,:) = 1.0/(this%lpul(i,:,:)*c**2) enddo this%cpul(size(this%segments)+1, :,:) = this%cpul(size(this%segments), :,:) From 6dcaa1f056e9983870a7448a61e8a98a644aed3c Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 10:26:17 +0100 Subject: [PATCH 03/10] Adds c_vacuum to parameters in fdetypes --- src_main_pub/fdetypes.F90 | 2 +- src_mtln/mtl.F90 | 9 ++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src_main_pub/fdetypes.F90 b/src_main_pub/fdetypes.F90 index 608de5e5..6e8f2fae 100755 --- a/src_main_pub/fdetypes.F90 +++ b/src_main_pub/fdetypes.F90 @@ -193,7 +193,7 @@ module FDETYPES 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12 real (KIND=RKIND), PARAMETER :: MU_VACUUM = & 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6 - + real (kind=rkind), parameter :: c_vacuum = 1.0_RKIND/sqrt(EPSILON_VACUUM*MU_VACUUM) REAL (KIND=RKIND_tiempo) :: dt0 !aqui para OLDrlo accesible en resume pscale diff --git a/src_mtln/mtl.F90 b/src_mtln/mtl.F90 index b4883e08..cae711e0 100644 --- a/src_mtln/mtl.F90 +++ b/src_mtln/mtl.F90 @@ -6,9 +6,9 @@ module mtl_mod use mtln_types_mod, only: segment_t, multipolar_expansion_t use multipolar_expansion_mod, only: getCellCapacitanceOnBox, getCellInductanceOnBox #ifdef CompileWithMPI - use fdetypes, only: SUBCOMM_MPI, REALSIZE, INTEGERSIZE, pi, mu_vacuum, epsilon_vacuum, RKIND_wires + use fdetypes, only: SUBCOMM_MPI, REALSIZE, INTEGERSIZE, pi, mu_vacuum, c_vacuum, RKIND_wires #else - use fdetypes, only: pi, mu_vacuum, epsilon_vacuum, RKIND_wires + use fdetypes, only: pi, mu_vacuum, c_vacuum, RKIND_wires #endif implicit none #ifdef CompileWithMPI @@ -269,11 +269,10 @@ subroutine computeLCParameters(this, multipolar_expansion) subroutine computeLCParametersFromRadius(this, rad) class(mtl_t) :: this real, intent(in) :: rad - REAL (KIND=RKIND_wires) :: invMu, c + REAL (KIND=RKIND_wires) :: invMu integer :: i real :: d1, d2 invMu = 1.0/mu_vacuum - c = 1/sqrt(epsilon_vacuum*mu_vacuum) do i = 1, size(this%segments) d1 = this%segments(i)%d1 d2 = this%segments(i)%d2 @@ -290,7 +289,7 @@ subroutine computeLCParametersFromRadius(this, rad) this%lpul(i,:,:) = this%lpul(i,:,:) & /(1.0_RKIND_wires-pi*rad**2.0_RKIND_wires /(d1*d2)) endif - this%cpul(i,:,:) = 1.0/(this%lpul(i,:,:)*c**2) + this%cpul(i,:,:) = 1.0/(this%lpul(i,:,:)*c_vacuum**2) enddo this%cpul(size(this%segments)+1, :,:) = this%cpul(size(this%segments), :,:) From 35caa1a6a2ff304ca28afaae56e9c14680a372a3 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 13:03:05 +0100 Subject: [PATCH 04/10] Holland test: Separates testing of wire without mtln, wire in mtln and unshielded multiwire in mtln --- test/pyWrapper/test_full_system.py | 65 ++++++++++++++++++------------ 1 file changed, 40 insertions(+), 25 deletions(-) diff --git a/test/pyWrapper/test_full_system.py b/test/pyWrapper/test_full_system.py index 5ecae5a3..8311bcea 100644 --- a/test/pyWrapper/test_full_system.py +++ b/test/pyWrapper/test_full_system.py @@ -174,48 +174,63 @@ def test_coated_antenna(tmp_path): p_solved['current_0'].to_numpy(), rtol=0.0, atol=10e-8) - -def test_holland(tmp_path): +@mtln_skip +def test_holland_single_wire(tmp_path): fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][0] = createWire(id = 1, r = 0.02) + solver.run() + p = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) + + expected_f = json.load(open(OUTPUTS_FOLDER+'holland1981.fdtd_mid_point_Wz_11_11_12_s2_12.json')) + expected_t, expected_i = np.array([]), np.array([]) + for data in expected_f['datasetColl'][0]['data']: + expected_t = np.append(expected_t, float(data['value'][0])) + expected_i = np.append(expected_i, float(data['value'][1])) - # if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - # solver['materials'][0] = createWire(id = 1, r = 0.02) - # elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - # solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) + expected_i_interp = np.interp(p['time']-3.05*1e-9, expected_t, expected_i) + assert np.allclose(expected_i_interp, -p['current'], rtol=1e-4, atol=5e-5) +@no_mtln_skip +def test_holland_wire(tmp_path): + fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' + solver = FDTD(input_filename=fn, + path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) solver['materials'][0] = createWire(id = 1, r = 0.02) solver.run() - p_wire = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) + p = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) + + expected_f = json.load(open(OUTPUTS_FOLDER+'holland1981.fdtd_mid_point_Wz_11_11_12_s2_12.json')) + expected_t, expected_i = np.array([]), np.array([]) + for data in expected_f['datasetColl'][0]['data']: + expected_t = np.append(expected_t, float(data['value'][0])) + expected_i = np.append(expected_i, float(data['value'][1])) + + expected_i_interp = np.interp(p['time']-3.05*1e-9, expected_t, expected_i) + assert np.allclose(expected_i_interp, p['current_0'], rtol=1e-4, atol=5e-5) +@no_mtln_skip +def test_holland_unshielded(tmp_path): + fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' + solver = FDTD(input_filename=fn, + path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) - solver.cleanUp() solver.run() - p_unshielded = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) - - + p = Probe(solver.getSolvedProbeFilenames("mid_point")[0]) + expected_f = json.load(open(OUTPUTS_FOLDER+'holland1981.fdtd_mid_point_Wz_11_11_12_s2_12.json')) expected_t, expected_i = np.array([]), np.array([]) for data in expected_f['datasetColl'][0]['data']: expected_t = np.append(expected_t, float(data['value'][0])) expected_i = np.append(expected_i, float(data['value'][1])) - expected_i_interp = np.interp(p_wire['time']-3.05*1e-9, expected_t, expected_i) - - # if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - # assert np.allclose( - # expected_i_interp, - # -p_solved['current'], - # rtol=1e-4, atol=5e-5) - # elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - # assert np.allclose( - # expected_i_interp, - # p_solved['current_0'], - # rtol=1e-4, atol=5e-5) - assert np.allclose(expected_i_interp, p_wire['current_0'], rtol=1e-4, atol=5e-5) - assert np.allclose(expected_i_interp, p_unshielded['current_0'], rtol=1e-4, atol=5e-5) + expected_i_interp = np.interp(p['time']-3.05*1e-9, expected_t, expected_i) + assert np.allclose(expected_i_interp, p['current_0'], rtol=1e-4, atol=5e-5) + @no_mtln_skip From 501dbd70bba5b54b2fb315b4cae71695b12689e9 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 13:21:40 +0100 Subject: [PATCH 05/10] Adds singlewire/wire/unshielded tests for towelHanger and lineIntegralProbe --- test/pyWrapper/test_full_system.py | 88 ++++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 17 deletions(-) diff --git a/test/pyWrapper/test_full_system.py b/test/pyWrapper/test_full_system.py index 8311bcea..c70642af 100644 --- a/test/pyWrapper/test_full_system.py +++ b/test/pyWrapper/test_full_system.py @@ -4,13 +4,35 @@ from sys import platform from scipy import signal -def test_lineIntegralProbe(tmp_path): +@mtln_skip +def test_lineIntegralProbe_single_wire(tmp_path): fn = CASES_FOLDER + 'lineIntegralProbe/lineIntegralProbe_plates.fdtd.json' solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][0] = createWire(id = 1, r = 0.001) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) + solver['materials'][0] = createWire(id = 1, r = 0.001) + solver.run() + + pf = 'lineIntegralProbe_plates.fdtd_vprobe_LI_20_20_10.dat' + li_probe = Probe(solver.getSolvedProbeFilenames("vprobe_LI_20_20_10")[0]) + expected = Probe(OUTPUTS_FOLDER+pf) + np.allclose(li_probe['lineIntegral'].to_numpy(), expected['lineIntegral'].to_numpy(), rtol =0.01 , atol=0.01) + +@no_mtln_skip +def test_lineIntegralProbe_wire(tmp_path): + fn = CASES_FOLDER + 'lineIntegralProbe/lineIntegralProbe_plates.fdtd.json' + solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][0] = createWire(id = 1, r = 0.001) + solver.run() + + pf = 'lineIntegralProbe_plates.fdtd_vprobe_LI_20_20_10.dat' + li_probe = Probe(solver.getSolvedProbeFilenames("vprobe_LI_20_20_10")[0]) + expected = Probe(OUTPUTS_FOLDER+pf) + np.allclose(li_probe['lineIntegral'].to_numpy(), expected['lineIntegral'].to_numpy(), rtol =0.01 , atol=0.01) + +@no_mtln_skip +def test_lineIntegralProbe_unshielded(tmp_path): + fn = CASES_FOLDER + 'lineIntegralProbe/lineIntegralProbe_plates.fdtd.json' + solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) solver.run() pf = 'lineIntegralProbe_plates.fdtd_vprobe_LI_20_20_10.dat' @@ -323,14 +345,31 @@ def test_unshielded_multiwires(tmp_path): p_solved['current_1'], rtol=1e-3, atol=0.01) +@mtln_skip def test_towelHanger(tmp_path): fn = CASES_FOLDER + 'towelHanger/towelHanger.fdtd.json' solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][0] = createWire(id = 1, r = 0.1e-3) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) + solver['materials'][0] = createWire(id = 1, r = 0.1e-3) + + p_solved = [Probe(solver.getSolvedProbeFilenames("wire_start")[0]), + Probe(solver.getSolvedProbeFilenames("wire_mid")[0]), + Probe(solver.getSolvedProbeFilenames("wire_end")[0])] + + p_expected = [Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_start_Wz_27_25_30_s1.dat'), + Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_mid_Wx_35_25_32_s5.dat'), + Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_end_Wz_43_25_30_s4.dat')] + + assert np.allclose(p_expected[0]['current'], p_solved[0]['current'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[1]['current'], p_solved[1]['current'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[2]['current'], p_solved[2]['current'], rtol=5e-3, atol=5e-4) + +@no_mtln_skip +def test_towelHanger(tmp_path): + fn = CASES_FOLDER + 'towelHanger/towelHanger.fdtd.json' + solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) + solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) solver.run() p_solved = [Probe(solver.getSolvedProbeFilenames("wire_start")[0]), @@ -341,14 +380,29 @@ def test_towelHanger(tmp_path): Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_mid_Wx_35_25_32_s5.dat'), Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_end_Wz_43_25_30_s4.dat')] - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - assert np.allclose(p_expected[0]['current'], p_solved[0]['current'], rtol=5e-3, atol=5e-4) - assert np.allclose(p_expected[1]['current'], p_solved[1]['current'], rtol=5e-3, atol=5e-4) - assert np.allclose(p_expected[2]['current'], p_solved[2]['current'], rtol=5e-3, atol=5e-4) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - assert np.allclose(p_expected[0]['current'], -p_solved[0]['current_0'], rtol=5e-3, atol=5e-4) - assert np.allclose(p_expected[1]['current'], -p_solved[1]['current_0'], rtol=5e-3, atol=5e-4) - assert np.allclose(p_expected[2]['current'], p_solved[2]['current_0'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[0]['current'], -p_solved[0]['current_0'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[1]['current'], -p_solved[1]['current_0'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[2]['current'], p_solved[2]['current_0'], rtol=5e-3, atol=5e-4) + +@no_mtln_skip +def test_towelHanger(tmp_path): + fn = CASES_FOLDER + 'towelHanger/towelHanger.fdtd.json' + solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) + solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) + solver.run() + + p_solved = [Probe(solver.getSolvedProbeFilenames("wire_start")[0]), + Probe(solver.getSolvedProbeFilenames("wire_mid")[0]), + Probe(solver.getSolvedProbeFilenames("wire_end")[0])] + + p_expected = [Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_start_Wz_27_25_30_s1.dat'), + Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_mid_Wx_35_25_32_s5.dat'), + Probe(OUTPUTS_FOLDER+'towelHanger.fdtd_wire_end_Wz_43_25_30_s4.dat')] + + assert np.allclose(p_expected[0]['current'], -p_solved[0]['current_0'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[1]['current'], -p_solved[1]['current_0'], rtol=5e-3, atol=5e-4) + assert np.allclose(p_expected[2]['current'], p_solved[2]['current_0'], rtol=5e-3, atol=5e-4) @no_hdf_skip From 9df33cd220ea9036bcb959840132533abdeed74b Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 13:36:35 +0100 Subject: [PATCH 06/10] Adds remaining tests for no mtln and mtln wire cases --- test/pyWrapper/test_full_system.py | 369 ++++++++++++++++++++++++++--- test/pyWrapper/test_integration.py | 48 +++- 2 files changed, 372 insertions(+), 45 deletions(-) diff --git a/test/pyWrapper/test_full_system.py b/test/pyWrapper/test_full_system.py index c70642af..e86820ea 100644 --- a/test/pyWrapper/test_full_system.py +++ b/test/pyWrapper/test_full_system.py @@ -534,39 +534,88 @@ def test_sgbc_shielding_effectiveness(tmp_path): assert np.allclose(fdtd_s21_db, anal_s21_db, rtol=0.05) -def test_sgbc_structured_resistance(tmp_path): +@mtln_skip +def test_sgbc_structured_resistance_single_wire(tmp_path): fn = CASES_FOLDER + 'sgbcResistance/sgbcResistance.fdtd.json' solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][2] = createWire(id = 3, r = 1e-4) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][2] = createUnshieldedWire(id = 3, lpul = 5.497210529384488e-07, cpul = 2.0212271390586895e-11) + solver['materials'][2] = createWire(id = 3, r = 1e-4) + solver.run() + + i = Probe(solver.getSolvedProbeFilenames("Bulk_probe")[0]).data['current'] + assert np.allclose(i.array[-101:-1], np.ones(100)*i.array[-100], rtol=1e-3) + assert np.allclose(-1/i.array[-101:-1], np.ones(100)*(50+45), rtol=0.05) +@no_mtln_skip +def test_sgbc_structured_resistance_wire(tmp_path): + fn = CASES_FOLDER + 'sgbcResistance/sgbcResistance.fdtd.json' + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + + solver['materials'][2] = createWire(id = 3, r = 1e-4) solver.run() i = Probe(solver.getSolvedProbeFilenames("Bulk_probe")[0]).data['current'] assert np.allclose(i.array[-101:-1], np.ones(100)*i.array[-100], rtol=1e-3) - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - assert np.allclose(-1/i.array[-101:-1], np.ones(100)*(50+45), rtol=0.05) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - assert np.allclose(1/i.array[-101:-1], np.ones(100)*(50+45), rtol=0.05) + assert np.allclose(1/i.array[-101:-1], np.ones(100)*(50+45), rtol=0.05) +@no_mtln_skip +def test_sgbc_structured_resistance_unshielded(tmp_path): + fn = CASES_FOLDER + 'sgbcResistance/sgbcResistance.fdtd.json' + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + + solver['materials'][2] = createUnshieldedWire(id = 3, lpul = 5.497210529384488e-07, cpul = 2.0212271390586895e-11) + solver.run() -def test_pec_overlapping_sgbcs(tmp_path): + i = Probe(solver.getSolvedProbeFilenames("Bulk_probe")[0]).data['current'] + assert np.allclose(i.array[-101:-1], np.ones(100)*i.array[-100], rtol=1e-3) + assert np.allclose(1/i.array[-101:-1], np.ones(100)*(50+45), rtol=0.05) + +@mtln_skip +def test_pec_overlapping_sgbcs_single_wire(tmp_path): """ Test that PEC surfaces overlapping SGBC surfaces prioritize PEC. """ fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json' # Runs case without overlap. solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][3] = createWire(id = 3, r = 1e-4) + solver.run() - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][3] = createWire(id = 3, r = 1e-4) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][3] = createUnshieldedWire(id = 3, lpul = 5.497210529384488e-07, cpul = 2.0212271390586895e-11) + p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0]) + t = p['time'].to_numpy() + iSGBC = p['current'].to_numpy() + # Adds current SGBC elements as PEC. Now both are defined over same surface. + sgbcElementIds = solver["materialAssociations"][1]["elementIds"] + solver['materialAssociations'][0]["elementIds"].extend(sgbcElementIds) + solver.cleanUp() solver.run() + iPEC = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])['current'].to_numpy() + + + # For debugging only. + # plt.figure() + # plt.plot(t, iSGBC,'.-', label='SGBC case') + # plt.plot(t, iPEC,'.-', label='PEC overlapping') + # plt.grid(which='both') + # plt.legend() + # plt.show() + + + # Checks values are different due to PEC prioritization. + assert np.all(np.greater(np.abs(iPEC[1000:]), np.abs(iSGBC[1000:]))) + +@no_mtln_skip +def test_pec_overlapping_sgbcs_wire(tmp_path): + """ Test that PEC surfaces overlapping SGBC surfaces prioritize PEC. + """ + fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json' + + # Runs case without overlap. + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][3] = createWire(id = 3, r = 1e-4) + solver.run() + p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0]) t = p['time'].to_numpy() iSGBC = p['current'].to_numpy() @@ -590,8 +639,43 @@ def test_pec_overlapping_sgbcs(tmp_path): # Checks values are different due to PEC prioritization. assert np.all(np.greater(np.abs(iPEC[1000:]), np.abs(iSGBC[1000:]))) + +@no_mtln_skip +def test_pec_overlapping_sgbcs_unshielded(tmp_path): + """ Test that PEC surfaces overlapping SGBC surfaces prioritize PEC. + """ + fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json' + + # Runs case without overlap. + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][3] = createUnshieldedWire(id = 3, lpul = 5.497210529384488e-07, cpul = 2.0212271390586895e-11) + solver.run() + + p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0]) + t = p['time'].to_numpy() + iSGBC = p['current'].to_numpy() + + # Adds current SGBC elements as PEC. Now both are defined over same surface. + sgbcElementIds = solver["materialAssociations"][1]["elementIds"] + solver['materialAssociations'][0]["elementIds"].extend(sgbcElementIds) + solver.cleanUp() + solver.run() + iPEC = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])['current'].to_numpy() + + + # For debugging only. + # plt.figure() + # plt.plot(t, iSGBC,'.-', label='SGBC case') + # plt.plot(t, iPEC,'.-', label='PEC overlapping') + # plt.grid(which='both') + # plt.legend() + # plt.show() -def test_sgbc_overlapping_sgbc(tmp_path): + + # Checks values are different due to PEC prioritization. + assert np.all(np.greater(np.abs(iPEC[1000:]), np.abs(iSGBC[1000:]))) +@mtln_skip +def test_sgbc_overlapping_sgbc_single_wire(tmp_path): """ Test that SGBC surfaces overlapping SGBC surfaces prioritize first in MatAss. """ fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json' @@ -599,11 +683,81 @@ def test_sgbc_overlapping_sgbc(tmp_path): # Runs case without overlap. solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) # Changes materialId in first SGBC in MatAss to material with larger conductivity. - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][3] = createWire(id = 3, r = 1e-4) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][3] = createUnshieldedWire(id = 3, lpul = 5.497210529384488e-07, cpul = 2.0212271390586895e-11) + solver['materials'][3] = createWire(id = 3, r = 1e-4) + solver['materialAssociations'][1]["materialId"] = 6 + solver.cleanUp() + solver.run() + p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0]) + + t = p['time'].to_numpy() + iSGBC_top = p['current'].to_numpy() + + # Changes materialId in second SGBC in MatAss to material with larger conductivity. + solver['materialAssociations'][1]["materialId"] = 2 + solver['materialAssociations'][2]["materialId"] = 6 + solver.cleanUp() + solver.run() + iSGBC_bottom = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])['current'].to_numpy() + + # For debugging only. + # plt.figure() + # plt.plot(t, iSGBC_top,'.-', label='SGBC sigma = 40 S/m, top') + # plt.plot(t, iSGBC_bottom,'.-', label='SGBC sigma = 20 S/m, bottom') + # plt.grid(which='both') + # plt.legend() + # plt.show() + + # Checks values are different due to prioritization of first written. + assert np.all(np.greater(np.abs(iSGBC_top[1000:]), np.abs(iSGBC_bottom[1000:]))) + +@no_mtln_skip +def test_sgbc_overlapping_sgbc_wire(tmp_path): + """ Test that SGBC surfaces overlapping SGBC surfaces prioritize first in MatAss. + """ + fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json' + # Runs case without overlap. + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + # Changes materialId in first SGBC in MatAss to material with larger conductivity. + solver['materials'][3] = createWire(id = 3, r = 1e-4) + solver['materialAssociations'][1]["materialId"] = 6 + solver.cleanUp() + solver.run() + p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0]) + + t = p['time'].to_numpy() + iSGBC_top = p['current'].to_numpy() + + # Changes materialId in second SGBC in MatAss to material with larger conductivity. + solver['materialAssociations'][1]["materialId"] = 2 + solver['materialAssociations'][2]["materialId"] = 6 + solver.cleanUp() + solver.run() + iSGBC_bottom = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])['current'].to_numpy() + + + # For debugging only. + # plt.figure() + # plt.plot(t, iSGBC_top,'.-', label='SGBC sigma = 40 S/m, top') + # plt.plot(t, iSGBC_bottom,'.-', label='SGBC sigma = 20 S/m, bottom') + # plt.grid(which='both') + # plt.legend() + # plt.show() + + + # Checks values are different due to prioritization of first written. + assert np.all(np.greater(np.abs(iSGBC_top[1000:]), np.abs(iSGBC_bottom[1000:]))) + +@no_mtln_skip +def test_sgbc_overlapping_sgbc_unshielded(tmp_path): + """ Test that SGBC surfaces overlapping SGBC surfaces prioritize first in MatAss. + """ + fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json' + + # Runs case without overlap. + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + # Changes materialId in first SGBC in MatAss to material with larger conductivity. + solver['materials'][3] = createUnshieldedWire(id = 3, lpul = 5.497210529384488e-07, cpul = 2.0212271390586895e-11) solver['materialAssociations'][1]["materialId"] = 6 solver.cleanUp() solver.run() @@ -745,17 +899,74 @@ def testCanExecuteFDTDFromFolderWithSpacesAndCanProcessAdditionalArguments(tmp_p assert (Probe(solver.getSolvedProbeFilenames("outside")[0]) is not None) assert (solver.getVTKMap()[0] is not None) - -def test_nodal_source(tmp_path): +@mtln_skip +def test_nodal_source_single_wire(tmp_path): fn = CASES_FOLDER + "nodalSource/nodalSource.fdtd.json" assert (os.path.isfile(fn)) solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][1] = createWire(id = 2, r = 0.1e-5, rpul=10000.0) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][1] = createUnshieldedWire(id = 2, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11, rpul = 10000.0) + solver['materials'][1] = createUnshieldedWire(id = 2, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11, rpul = 10000.0) + solver.run() + resistanceBulkProbe = Probe( \ + solver.getSolvedProbeFilenames("Bulk probe Resistance")[0]) + nodalBulkProbe = Probe( \ + solver.getSolvedProbeFilenames("Bulk probe Nodal Source")[0]) + excitation = ExcitationFile( \ + excitation_filename=solver.getExcitationFile("predefinedExcitation")[0]) + + # For debugging. + # plt.figure() + # plt.plot(resistanceBulkProbe['time'].to_numpy(), + # resistanceBulkProbe['current'].to_numpy(), label='BP Current@resistance') + # plt.plot(excitation.data['time'].to_numpy(), + # excitation.data['value'].to_numpy(), label='excited current') + # plt.plot(nodalBulkProbe['time'].to_numpy(), + # -nodalBulkProbe['current'].to_numpy(), label='BP Current@nodal source') + # plt.legend() + + exc = np.interp(nodalBulkProbe['time'].to_numpy(), + excitation.data['time'].to_numpy(), + excitation.data['value'].to_numpy()) + assert np.corrcoef(exc, -nodalBulkProbe['current'])[0,1] > 0.999 + assert np.corrcoef(-nodalBulkProbe['current'], resistanceBulkProbe['current'])[0,1] > 0.998 + +@no_mtln_skip +def test_nodal_source_wire(tmp_path): + fn = CASES_FOLDER + "nodalSource/nodalSource.fdtd.json" + assert (os.path.isfile(fn)) + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][1] = createWire(id = 2, r = 0.1e-5, rpul=10000.0) + solver.run() + resistanceBulkProbe = Probe( \ + solver.getSolvedProbeFilenames("Bulk probe Resistance")[0]) + nodalBulkProbe = Probe( \ + solver.getSolvedProbeFilenames("Bulk probe Nodal Source")[0]) + excitation = ExcitationFile( \ + excitation_filename=solver.getExcitationFile("predefinedExcitation")[0]) + + # For debugging. + # plt.figure() + # plt.plot(resistanceBulkProbe['time'].to_numpy(), + # resistanceBulkProbe['current'].to_numpy(), label='BP Current@resistance') + # plt.plot(excitation.data['time'].to_numpy(), + # excitation.data['value'].to_numpy(), label='excited current') + # plt.plot(nodalBulkProbe['time'].to_numpy(), + # -nodalBulkProbe['current'].to_numpy(), label='BP Current@nodal source') + # plt.legend() + + exc = np.interp(nodalBulkProbe['time'].to_numpy(), + excitation.data['time'].to_numpy(), + excitation.data['value'].to_numpy()) + assert np.corrcoef(exc, -nodalBulkProbe['current'])[0,1] > 0.999 + assert np.corrcoef(-nodalBulkProbe['current'], resistanceBulkProbe['current'])[0,1] > 0.998 + +@no_mtln_skip +def test_nodal_source_unshielded(tmp_path): + fn = CASES_FOLDER + "nodalSource/nodalSource.fdtd.json" + assert (os.path.isfile(fn)) + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + solver['materials'][1] = createUnshieldedWire(id = 2, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11, rpul = 10000.0) solver.run() resistanceBulkProbe = Probe( \ @@ -1196,18 +1407,14 @@ def test_negative_offset_in_x(tmp_path): assert np.corrcoef(probeL['current'].to_numpy(), I_interp)[0, 1] > 0.999 assert np.allclose(probeR['current'].to_numpy(), 0.0, atol=3e-3) -@mtln_skip -def test_conformal_impedance_cylinder(tmp_path): +@mtln_skip +def test_conformal_impedance_cylinder_single_wire(tmp_path): case_name = 'conformal_impedance_cylinder_conformal' solver = FDTD(input_filename=TEST_DATA_FOLDER+'cases/conformal_impedance_cylinder/'+case_name+'.fdtd.json', path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) solver.cleanUp() - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][2] = createWire(id = 3, r = 0.1e-3) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][2] = createUnshieldedWire(id = 3, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) - + solver['materials'][2] = createWire(id = 3, r = 0.1e-3) solver.run() assert solver.hasFinishedSuccessfully() bulk_conf = Probe(solver.getSolvedProbeFilenames("BulkProbe")[0]) @@ -1217,11 +1424,101 @@ def test_conformal_impedance_cylinder(tmp_path): run_in_folder=tmp_path) solver.cleanUp() - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][2] = createWire(id = 3, r = 0.1e-3) - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][2] = createUnshieldedWire(id = 3, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) + solver['materials'][2] = createWire(id = 3, r = 0.1e-3) + solver.run() + assert solver.hasFinishedSuccessfully() + bulk = Probe(solver.getSolvedProbeFilenames("BulkProbe")[0]) + + #discrete fourier transforms + exc_file = "predefinedExcitation.1.exc" + exc = pd.read_csv(exc_file, sep='\\s+') + exc = exc.rename(columns={ + exc.columns[0]: 'time', + exc.columns[1]: 'V' + }) + new_freqs = np.geomspace(1e3, 1e7, num=100) + Vexc = exc["V"].to_numpy() + texc = exc["time"].to_numpy() + dt_exc = texc[1]-texc[0] + Vfexc = dt_exc*np.array([np.sum(Vexc * np.exp(-1j * 2 * np.pi * f * texc)) for f in new_freqs]) + + Ibulk = bulk["current"].to_numpy() + tbulk = bulk["time"].to_numpy() + dt_bulk = tbulk[1]-tbulk[0] + Ifbulk = dt_bulk*np.array([np.sum(Ibulk * np.exp(-1j * 2 * np.pi * f * tbulk)) for f in new_freqs]) + + Ibulk_conf = bulk_conf["current"].to_numpy() + tbulk_conf = bulk_conf["time"].to_numpy() + dt_bulk_conf = tbulk_conf[1]-tbulk_conf[0] + Ifbulk_conf = dt_bulk_conf*np.array([np.sum(Ibulk_conf * np.exp(-1j * 2 * np.pi * f * tbulk_conf)) for f in new_freqs]) + + #impedance comparison + assert np.allclose(np.abs(Vfexc/Ifbulk), np.abs(Vfexc/Ifbulk_conf), rtol=0.01, atol=0.001) + +@no_mtln_skip +def test_conformal_impedance_cylinder_wire(tmp_path): + case_name = 'conformal_impedance_cylinder_conformal' + solver = FDTD(input_filename=TEST_DATA_FOLDER+'cases/conformal_impedance_cylinder/'+case_name+'.fdtd.json', path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) + solver.cleanUp() + + solver['materials'][2] = createWire(id = 3, r = 0.1e-3) + solver.run() + assert solver.hasFinishedSuccessfully() + bulk_conf = Probe(solver.getSolvedProbeFilenames("BulkProbe")[0]) + + case_name = 'conformal_impedance_cylinder_staircase' + solver = FDTD(input_filename=TEST_DATA_FOLDER+'cases/conformal_impedance_cylinder/'+case_name+'.fdtd.json', path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) + solver.cleanUp() + + solver['materials'][2] = createWire(id = 3, r = 0.1e-3) + solver.run() + assert solver.hasFinishedSuccessfully() + bulk = Probe(solver.getSolvedProbeFilenames("BulkProbe")[0]) + + #discrete fourier transforms + exc_file = "predefinedExcitation.1.exc" + exc = pd.read_csv(exc_file, sep='\\s+') + exc = exc.rename(columns={ + exc.columns[0]: 'time', + exc.columns[1]: 'V' + }) + new_freqs = np.geomspace(1e3, 1e7, num=100) + Vexc = exc["V"].to_numpy() + texc = exc["time"].to_numpy() + dt_exc = texc[1]-texc[0] + Vfexc = dt_exc*np.array([np.sum(Vexc * np.exp(-1j * 2 * np.pi * f * texc)) for f in new_freqs]) + Ibulk = bulk["current"].to_numpy() + tbulk = bulk["time"].to_numpy() + dt_bulk = tbulk[1]-tbulk[0] + Ifbulk = dt_bulk*np.array([np.sum(Ibulk * np.exp(-1j * 2 * np.pi * f * tbulk)) for f in new_freqs]) + + Ibulk_conf = bulk_conf["current"].to_numpy() + tbulk_conf = bulk_conf["time"].to_numpy() + dt_bulk_conf = tbulk_conf[1]-tbulk_conf[0] + Ifbulk_conf = dt_bulk_conf*np.array([np.sum(Ibulk_conf * np.exp(-1j * 2 * np.pi * f * tbulk_conf)) for f in new_freqs]) + + #impedance comparison + assert np.allclose(np.abs(Vfexc/Ifbulk), np.abs(Vfexc/Ifbulk_conf), rtol=0.01, atol=0.001) + +@no_mtln_skip +def test_conformal_impedance_cylinder(tmp_path): + case_name = 'conformal_impedance_cylinder_conformal' + solver = FDTD(input_filename=TEST_DATA_FOLDER+'cases/conformal_impedance_cylinder/'+case_name+'.fdtd.json', path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) + solver.cleanUp() + solver['materials'][2] = createUnshieldedWire(id = 3, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) + solver.run() + assert solver.hasFinishedSuccessfully() + bulk_conf = Probe(solver.getSolvedProbeFilenames("BulkProbe")[0]) + + case_name = 'conformal_impedance_cylinder_staircase' + solver = FDTD(input_filename=TEST_DATA_FOLDER+'cases/conformal_impedance_cylinder/'+case_name+'.fdtd.json', path_to_exe=SEMBA_EXE, + run_in_folder=tmp_path) + solver.cleanUp() + solver['materials'][2] = createUnshieldedWire(id = 3, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11) solver.run() assert solver.hasFinishedSuccessfully() bulk = Probe(solver.getSolvedProbeFilenames("BulkProbe")[0]) diff --git a/test/pyWrapper/test_integration.py b/test/pyWrapper/test_integration.py index 3ffab0eb..fadeaf6e 100644 --- a/test/pyWrapper/test_integration.py +++ b/test/pyWrapper/test_integration.py @@ -1,20 +1,50 @@ from utils import * -def test_holland_case_checking_number_of_outputs(tmp_path): +@mtln_skip +def test_holland_case_checking_number_of_outputs_single_wire(tmp_path): fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) number_of_steps = 10 solver['general']['numberOfSteps'] = number_of_steps - if (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "OFF"): - solver['materials'][0] = createWire(id = 1, r = 0.02) - outfile = 'holland1981.fdtd_mid_point_Wz_11_11_12_s2.dat' - elif (os.getenv("SEMBA_FDTD_ENABLE_MTLN") == "ON"): - solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) - outfile = 'holland1981.fdtd_mid_point_single_wire_I_11_11_12.dat' - - + solver['materials'][0] = createWire(id = 1, r = 0.02) + outfile = 'holland1981.fdtd_mid_point_Wz_11_11_12_s2.dat' + solver.run() + + probe_files = solver.getSolvedProbeFilenames("mid_point") + + assert len(probe_files) == 1 + assert outfile == probe_files[0] + assert countLinesInFile(probe_files[0]) == number_of_steps + 2 + +@no_mtln_skip +def test_holland_case_checking_number_of_outputs_wire(tmp_path): + fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + + number_of_steps = 10 + solver['general']['numberOfSteps'] = number_of_steps + + solver['materials'][0] = createWire(id = 1, r = 0.02) + outfile = 'holland1981.fdtd_mid_point_single_wire_I_11_11_12.dat' + solver.run() + + probe_files = solver.getSolvedProbeFilenames("mid_point") + + assert len(probe_files) == 1 + assert outfile == probe_files[0] + assert countLinesInFile(probe_files[0]) == number_of_steps + 2 + +@no_mtln_skip +def test_holland_case_checking_number_of_outputs_unshielded(tmp_path): + fn = CASES_FOLDER + 'holland/holland1981.fdtd.json' + solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) + + number_of_steps = 10 + solver['general']['numberOfSteps'] = number_of_steps + solver['materials'][0] = createUnshieldedWire(id = 1, lpul = 6.52188703e-08, cpul = 1.7060247700000001e-10) + outfile = 'holland1981.fdtd_mid_point_single_wire_I_11_11_12.dat' solver.run() probe_files = solver.getSolvedProbeFilenames("mid_point") From 053ae6d94c0dd491579c8baa326fc75cef8f2061 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 13:50:16 +0100 Subject: [PATCH 07/10] removes change commited by error --- test/pyWrapper/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pyWrapper/utils.py b/test/pyWrapper/utils.py index b86eb333..7f1d1857 100644 --- a/test/pyWrapper/utils.py +++ b/test/pyWrapper/utils.py @@ -34,7 +34,7 @@ # Use of absolute path to avoid conflicts when changing directory. if platform == "linux": - SEMBA_EXE = os.path.join(os.getcwd(), 'build-dbg', 'bin', 'semba-fdtd') + SEMBA_EXE = os.path.join(os.getcwd(), 'build', 'bin', 'semba-fdtd') elif platform == "win32": SEMBA_EXE = os.path.join(os.getcwd(), 'build', 'bin', 'semba-fdtd.exe') From b6ef3da7d943957901e1a04b3f641ed1d0604b22 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 15:53:03 +0100 Subject: [PATCH 08/10] Minor the read rpul and gpul in different formats from wires and multiwires --- src_json_parser/smbjson.F90 | 33 ++++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/src_json_parser/smbjson.F90 b/src_json_parser/smbjson.F90 index b067cac1..b076cefc 100644 --- a/src_json_parser/smbjson.F90 +++ b/src_json_parser/smbjson.F90 @@ -71,6 +71,7 @@ module smbjson procedure, private :: getMatrixAt procedure, private :: getStrAt procedure, private :: existsAt + procedure, private :: dimensionAt procedure, private :: getDomain procedure, private :: buildPECPMCRegions procedure, private :: getMaterialAssociations @@ -2526,7 +2527,6 @@ function readMTLN(this) result (mtln_res) type(materialAssociation_t), dimension(:), allocatable :: cables class(cable_t), pointer :: ptr, read_cable integer :: i - cables = this%getMaterialAssociations([ & J_MAT_TYPE_SHIELDED_MULTIWIRE//' ',& J_MAT_TYPE_UNSHIELDED_MULTIWIRE ,& @@ -2534,7 +2534,6 @@ function readMTLN(this) result (mtln_res) ! spaces are needed to make strings have same length. ! Why? Because of FORTRAN! It only accepts fixed length strings for arrays. - mtln_res%connectors => readConnectors() call addConnIdToConnectorMap(connIdToConnector, mtln_res%connectors) if (size(cables) == 0) then @@ -3496,7 +3495,6 @@ function readMTLNCable(j_cable, despl) result(res) case default call WarnErrReport("Error reading cable: material type is not valid", .true.) end select - res%initial_connector => findConnectorWithId(j_cable%initialConnectorId) res%end_connector => findConnectorWithId(j_cable%endConnectorId) res%name = j_cable%name @@ -3583,6 +3581,7 @@ subroutine assignInCellProperties(res, mat, n) "Unshielded multiwires in cell properties must be defined by fixed OR multipolarExpansions, but not both.", .true.) end if end if + if (areFixedInCell) then res%cell_inductance_per_meter = this%getMatrixAt(mat%p, J_MAT_MULTIWIRE_INDUCTANCE,found) res%cell_capacitance_per_meter = this%getMatrixAt(mat%p, J_MAT_MULTIWIRE_CAPACITANCE,found) @@ -3602,23 +3601,27 @@ subroutine assignInCellProperties(res, mat, n) end if if (this%existsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE)) then - r = this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found) - if (found) then - res%resistance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found)) + m = this%dimensionAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE) + if (m == 0) then + allocate(r(1)) + r(1) = this%getRealAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found) else - res%resistance_per_meter = null_matrix + r = this%getRealsAt(mat%p, J_MAT_MULTIWIRE_RESISTANCE,found) end if + res%resistance_per_meter = vectorToDiagonalMatrix(r) else res%resistance_per_meter = null_matrix end if if (this%existsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE)) then - c = this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found) - if (found) then - res%conductance_per_meter = vectorToDiagonalMatrix(this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found)) + m = this%dimensionAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE) + if (m == 0) then + allocate(c(1)) + c(1) = this%getRealAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found) else - res%conductance_per_meter = null_matrix + c = this%getRealsAt(mat%p, J_MAT_MULTIWIRE_CONDUCTANCE,found) end if + res%conductance_per_meter = vectorToDiagonalMatrix(c) else res%conductance_per_meter = null_matrix end if @@ -4054,6 +4057,14 @@ function existsAt(this, place, path) result(res) call this%core%info(place, path, found=res) end function + function dimensionAt(this, place, path) result(res) + integer :: res + class(parser_t) :: this + type(json_value), pointer :: place + character(len=*) :: path + call this%core%info(place, path, n_children=res) + end function + function jsonValueFilterByKeyValues(this, srcs, key, values) result (res) class(parser_t) :: this type(json_value_ptr), dimension(:), allocatable :: res From 00206cba6612d0a34a699394fe63021005db593b Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 19 Feb 2026 16:09:52 +0100 Subject: [PATCH 09/10] Fixes error in material assignment in test --- test/pyWrapper/test_full_system.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pyWrapper/test_full_system.py b/test/pyWrapper/test_full_system.py index e86820ea..4288b065 100644 --- a/test/pyWrapper/test_full_system.py +++ b/test/pyWrapper/test_full_system.py @@ -904,7 +904,7 @@ def test_nodal_source_single_wire(tmp_path): fn = CASES_FOLDER + "nodalSource/nodalSource.fdtd.json" assert (os.path.isfile(fn)) solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path) - solver['materials'][1] = createUnshieldedWire(id = 2, lpul = 6.5183032590978384e-07, cpul = 1.7046017451862063e-11, rpul = 10000.0) + solver['materials'][1] = createWire(id = 2, r = 0.1e-5, rpul=10000.0) solver.run() resistanceBulkProbe = Probe( \ From 0aeb5e1836052c5ae15638ceb6f9092044eec356 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Fri, 20 Feb 2026 13:48:21 +0100 Subject: [PATCH 10/10] Now preallocates when reading sources for circuits. --- src_mtln/circuit.F90 | 45 +++++++++++++++++++++++++++++++++----------- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/src_mtln/circuit.F90 b/src_mtln/circuit.F90 index c11bd1ae..6cf40fe2 100644 --- a/src_mtln/circuit.F90 +++ b/src_mtln/circuit.F90 @@ -134,18 +134,41 @@ subroutine init(this, names, sources, netlist) type(source_t) function setSource(source_path) result(res) character(*), intent(in) :: source_path real :: time, value - integer :: io - allocate(res%time(0), res%value(0)) - if (source_path /= "" ) then - res%has_source = .true. - open(unit = 1, file = source_path) - do - read(1, *, iostat = io) time, value - if (io /= 0) exit - res%time = [res%time, time] - res%value = [res%value, value] - end do + integer :: io, line_count, i + + if (source_path == "" ) then + allocate(res%time(0), res%value(0)) + res%has_source = .false. + return end if + + res%has_source = .true. + + ! First pass: count the number of lines + line_count = 0 + open(unit = 1, file = source_path) + do + read(1, *, iostat = io) time, value + if (io /= 0) exit + line_count = line_count + 1 + end do + close(1) + + ! Allocate arrays with the exact size needed + allocate(res%time(line_count)) + allocate(res%value(line_count)) + + ! Second pass: fill the arrays + open(unit = 1, file = source_path) + do i = 1, line_count + read(1, *, iostat = io) res%time(i), res%value(i) + if (io /= 0) then + ! Handle unexpected read error + close(1) + error stop "Error reading excitation file" + end if + end do + close(1) end function subroutine loadNetlist(this, netlist)