Skip to content

260 parse wires when compiled with mtln #261

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Alberto-o merged 10 commits into from
Feb 20, 2026
Merged

260 parse wires when compiled with mtln #261

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
02f43d1
WIP
Alberto-o Feb 18, 2026
00c219b
Bouyated inductance model implemented for wires in mtln
Alberto-o Feb 19, 2026
6dcaa1f
Adds c_vacuum to parameters in fdetypes
Alberto-o Feb 19, 2026
35caa1a
Holland test: Separates testing of wire without mtln, wire in mtln an…
Alberto-o Feb 19, 2026
501dbd7
Adds singlewire/wire/unshielded tests for towelHanger and lineIntegra…
Alberto-o Feb 19, 2026
9df33cd
Adds remaining tests for no mtln and mtln wire cases
Alberto-o Feb 19, 2026
053ae6d
removes change commited by error
Alberto-o Feb 19, 2026
b6ef3da
Minor the read rpul and gpul in different formats from wires and mult…
Alberto-o Feb 19, 2026
00206cb
Fixes error in material assignment in test
Alberto-o Feb 19, 2026
0aeb5e1
Now preallocates when reading sources for circuits.
lmdiazangulo Feb 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src_json_parser/parser_tools.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,6 +191,7 @@ function getPixelFromElementId(mesh, id) result(res)

end function


function vectorToDiagonalMatrix(vector) result(res)
real, dimension(:), intent(in) :: vector
real, dimension(:, :), allocatable :: res
Expand Down
101 changes: 76 additions & 25 deletions src_json_parser/smbjson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -1225,8 +1226,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
Expand Down Expand Up @@ -2518,22 +2524,16 @@ 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.


mtln_res%connectors => readConnectors()
call addConnIdToConnectorMap(connIdToConnector, mtln_res%connectors)
if (size(cables) == 0) then
Expand Down Expand Up @@ -2589,6 +2589,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
Expand All @@ -2612,6 +2614,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
Expand Down Expand Up @@ -2712,7 +2716,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)
Expand Down Expand Up @@ -3264,7 +3269,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))
Expand Down Expand Up @@ -3303,7 +3309,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))
Expand Down Expand Up @@ -3336,7 +3343,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))
Expand Down Expand Up @@ -3478,7 +3486,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)
Expand All @@ -3487,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
Expand Down Expand Up @@ -3548,46 +3555,73 @@ 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
real, dimension(:), allocatable :: r, c

allocate(null_matrix(n,n), source = 0.0)

areFixedInCell = &
this%existsAt(mat%p, J_MAT_MULTIWIRE_INDUCTANCE) .and. &
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
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
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
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
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
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
Expand Down Expand Up @@ -3679,14 +3713,23 @@ 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)
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)
Expand Down Expand Up @@ -4014,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
Expand Down
2 changes: 1 addition & 1 deletion src_main_pub/fdetypes.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
45 changes: 34 additions & 11 deletions src_mtln/circuit.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
42 changes: 38 additions & 4 deletions src_mtln/mtl.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, c_vacuum, RKIND_wires
#else
use fdetypes, only: pi, mu_vacuum, c_vacuum, RKIND_wires
#endif
implicit none
#ifdef CompileWithMPI
Expand Down Expand Up @@ -59,6 +61,7 @@ module mtl_mod
procedure :: checkTimeStep
procedure :: allocatePULMatrices
procedure :: computeLCParameters
procedure :: computeLCParametersFromRadius
procedure :: initLC
procedure :: initRG
procedure :: initDirections
Expand Down Expand Up @@ -163,15 +166,15 @@ 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
character(len=*), intent(in) :: name
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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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,:,:)*c_vacuum**2)
enddo
this%cpul(size(this%segments)+1, :,:) = this%cpul(size(this%segments), :,:)

end subroutine


subroutine initDirections(this)
class(mtl_t) :: this
Expand Down
Loading
Loading