Modifies cable_t segment_relative_positions to hold external fields. WIP

in wires_mtln.f90

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required (VERSION 3.18)
 
-ENABLE_LANGUAGE(Fortran)
+enable_language (CXX Fortran)
 project(semba-fdtd Fortran)
 
 set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
@@ -64,8 +64,9 @@ if (CompileExecutable)
 		"src_main_pub/version.F90"
 		"src_wires_pub/wires_types.F90"
 		"src_wires_pub/wires.F90"
+		"src_wires_pub/wires_mtln.F90"
 	)
-	target_link_libraries(semba-fdtd semba-types smbjson)
+	target_link_libraries(semba-fdtd semba-types smbjson mtlnsolver ngspice_interface ${NGSPICE_LIB} stdc++)
 	include_directories(${CMAKE_BINARY_DIR}/mod)
 endif()
 

src_json_parser/smbjson.F90

@@ -1912,7 +1912,7 @@ function readMTLNCable(j_cable, is_read) result(res)
          end if
 
          res%step_size = buildStepSize(j_cable)
-         res%segment_relative_positions = mapSegmentsToGridCoordinates(j_cable)
+         res%external_field_segments = mapSegmentsToGridCoordinates(j_cable)
          ! write(*,*) 'id: ', this%getIntAt(j_cable, J_MATERIAL_ID, found)
          material = this%matTable%getId(this%getIntAt(j_cable, J_MATERIAL_ID, found))
          if (.not. found) &
@@ -2075,7 +2075,7 @@ subroutine assignPULProperties(res, mat, n)
 
       function mapSegmentsToGridCoordinates(j_cable) result(res)
          type(json_value), pointer :: j_cable
-         type(segment_relative_position_t), dimension(:), allocatable :: res
+         type(external_field_segment_t), dimension(:), allocatable :: res
          integer, dimension(:), allocatable :: elemIds
          type(polyline_t) :: p_line
 
@@ -2104,14 +2104,16 @@ function mapSegmentsToGridCoordinates(j_cable) result(res)
 
       function mapNegativeSegment(c1, c2) result(res)
          type(coordinate_t), intent(in) :: c1, c2
-         type(segment_relative_position_t) :: curr_pos
+         type(external_field_segment_t) :: curr_pos
          integer :: axis, i, n_segments
-         type(segment_relative_position_t), dimension(:), allocatable :: res
+         type(external_field_segment_t), dimension(:), allocatable :: res
 
          axis = findDirection(c2-c1)
          n_segments = abs(ceiling(c2%position(axis)) - floor(c1%position(axis)))
          allocate(res(n_segments))
          curr_pos%position = [(c1%position(i), i = 1, 3)]
+         curr_pos%Efield_main2wire => null()
+         curr_pos%Efield_wire2main => null()
 
          res = [(curr_pos, i = 1, n_segments)]
          res(:)%position(axis) = [(res(i)%position(axis) - i, i = 1, n_segments)]
@@ -2120,15 +2122,17 @@ function mapNegativeSegment(c1, c2) result(res)
 
       function mapPositiveSegment(c1, c2) result(res)
          type(coordinate_t), intent(in) :: c1, c2
-         type(segment_relative_position_t) :: curr_pos
+         type(external_field_segment_t) :: curr_pos
          integer :: axis, orientation, i, n_segments
-         type(segment_relative_position_t), dimension(:), allocatable :: res
+         type(external_field_segment_t), dimension(:), allocatable :: res
 
          axis = findDirection(c2-c1)
 
          n_segments = abs(floor(c2%position(axis)) - ceiling(c1%position(axis)))
          allocate(res(n_segments))
          curr_pos%position = [(c1%position(i), i = 1, 3)]
+         curr_pos%Efield_main2wire => null()
+         curr_pos%Efield_wire2main => null()
 
          res = [(curr_pos, i = 1, n_segments)]
          res(:)%position(axis) = [(res(i)%position(axis) + (i-1), i = 1, n_segments)]

src_main_pub/timestepping.F90

@@ -1367,7 +1367,7 @@ subroutine launch_simulation(sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx
          endif
 #endif
 #ifdef CompileWithWires_mtln  
-         if (thereAre%MTLNbundles) call AdvanceWiresE_mtln(sgg,Ex,Ey,Ez,Idxe,Idye,Idze,Idxh,Idyh,Idzh,eps0,mu0,mtln_solver)  
+         if (thereAre%MTLNbundles) call AdvanceWiresE_mtln(sgg, Idxe,Idye,Idze,Idxh,Idyh,Idzh,eps0,mu0,mtln_solver)  
 #endif 
          If (Thereare%PMLbodies) then !waveport absorbers
             call AdvancePMLbodyE

src_mtln/mtl.F90

@@ -3,7 +3,7 @@ module mtl_mod
     ! use NFDETypes
     use utils_mod
     use dispersive_mod
-    use mtln_types_mod, only: segment_relative_position_t
+    use mtln_types_mod, only: external_field_segment_t
 
     implicit none
 
@@ -20,7 +20,7 @@ module mtl_mod
         integer :: conductor_in_parent
 
         type(transfer_impedance_per_meter_t) :: transfer_impedance
-        type(segment_relative_position_t), allocatable, dimension(:) :: segment_relative_positions
+        type(external_field_segment_t), allocatable, dimension(:) :: external_field_segments
 
     contains
         procedure :: setTimeStep
@@ -84,7 +84,7 @@ function mtlHomogeneous(lpul, cpul, rpul, gpul, &
                             step_size, name, &
                             dt, parent_name, conductor_in_parent, &
                             transfer_impedance, &
-                            segment_relative_positions) result(res)
+                            external_field_segments) result(res)
         type(mtl_t) :: res
         real, intent(in), dimension(:,:) :: lpul, cpul, rpul, gpul
         real, intent(in), dimension(:) :: step_size
@@ -95,7 +95,7 @@ function mtlHomogeneous(lpul, cpul, rpul, gpul, &
         character(len=*), intent(in), optional :: parent_name
         integer, intent(in), optional :: conductor_in_parent
         type(transfer_impedance_per_meter_t), intent(in), optional :: transfer_impedance
-        type(segment_relative_position_t), intent(in), dimension(:), optional :: segment_relative_positions
+        type(external_field_segment_t), intent(in), dimension(:), optional :: external_field_segments
 
         res%name = name
         res%step_size =  step_size
@@ -130,15 +130,16 @@ function mtlHomogeneous(lpul, cpul, rpul, gpul, &
             res%transfer_impedance = transfer_impedance
         end if
 
-        if (present(segment_relative_positions)) then 
-            res%segment_relative_positions = segment_relative_positions
+        if (present(external_field_segments)) then 
+            res%external_field_segments = external_field_segments
         end if
     end function
 
     function mtlInhomogeneous(lpul, cpul, rpul, gpul, &
                             step_size, name, &
                             dt, parent_name, conductor_in_parent, &
-                            transfer_impedance) result(res)
+                            transfer_impedance, &
+                            external_field_segments) result(res)
         type(mtl_t) :: res
         real, intent(in), dimension(:,:,:) :: lpul, cpul, rpul, gpul
         real, intent(in), dimension(:) :: step_size
@@ -148,6 +149,7 @@ function mtlInhomogeneous(lpul, cpul, rpul, gpul, &
         character(len=*), intent(in), optional :: parent_name
         integer, intent(in), optional :: conductor_in_parent
         type(transfer_impedance_per_meter_t), intent(in), optional :: transfer_impedance
+        type(external_field_segment_t), intent(in), dimension(:), optional :: external_field_segments
 
         res%name = name
         res%step_size = step_size
@@ -181,6 +183,10 @@ function mtlInhomogeneous(lpul, cpul, rpul, gpul, &
             res%transfer_impedance = transfer_impedance
         end if
 
+        if (present(external_field_segments)) then 
+            res%external_field_segments = external_field_segments
+        end if
+
     end function
 
     subroutine initDirections(this)

src_mtln/mtl_bundle.F90

@@ -21,7 +21,7 @@ module mtl_bundle_mod
         real, dimension(:,:,:), allocatable :: v_term, i_term
         real, dimension(:,:,:), allocatable :: v_diff, i_diff
 
-        type(segment_relative_position_t), dimension(:), allocatable :: segment_relative_positions
+        type(external_field_segment_t), dimension(:), allocatable :: external_field_segments
 
 
     contains
@@ -37,7 +37,6 @@ module mtl_bundle_mod
         procedure :: advanceCurrent => bundle_advanceCurrent
         procedure :: addTransferImpedance => bundle_addTransferImpedance
         ! procedure :: setConnectorTransferImpedance
-        procedure :: setExternalCurrent => bundle_setExternalCurrent
         procedure :: setExternalVoltage => bundle_setExternalVoltage
         procedure :: updateExternalCurrent => bundle_updateExternalCurrent
 
@@ -64,7 +63,7 @@ function mtldCtor(levels, name) result(res)
         res%dt = levels(1)%lines(1)%dt
         res%step_size = levels(1)%lines(1)%step_size
         res%number_of_divisions = size(res%step_size,1)
-        res%segment_relative_positions = levels(1)%lines(1)%segment_relative_positions
+        res%external_field_segments = levels(1)%lines(1)%external_field_segments
         call res%initialAllocation()
         call res%mergePULMatrices(levels)
         call res%mergeDispersiveMatrices(levels)
@@ -309,23 +308,12 @@ subroutine bundle_advanceCurrent(this)
         ! call this%transfer_impedance%updatePhi(i_prev, i_now)
     end subroutine
 
-    subroutine bundle_setExternalCurrent(this, current)
+    subroutine bundle_setExternalVoltage(this)
         class(mtl_bundle_t) :: this
-        real, dimension(:), intent(in) :: current
-        this%i(1,:) = current(:)
-    end subroutine
-
-    subroutine bundle_setExternalVoltage(this, voltage)
-        class(mtl_bundle_t) :: this
-        real, dimension(:,:,:), intent(in) :: voltage
-        ! real, dimension(:), intent(in) :: voltage
-        integer, dimension(:), allocatable :: position
         integer :: i
         do i = 1, size(this%v,2) 
-            position = this%segment_relative_positions(i)%position
-            this%v(1, i) = voltage(position(1), position(2), position(3))
+            this%v(1, i) = this%external_field_segments(i)%Efield_main2wire * this%step_size(i)
         end do
-        ! this%v(1,:) = voltage(:)
     end subroutine
 
     subroutine bundle_updateExternalCurrent(this, current)
@@ -335,7 +323,7 @@ subroutine bundle_updateExternalCurrent(this, current)
         integer, dimension(:), allocatable :: position
         integer :: i
         do i = 1, size(this%i,2)
-            position = this%segment_relative_positions(i)%position
+            position = this%external_field_segments(i)%position
             current(position(1), position(2), position(3)) = this%i(1,i)
         end do
         ! current(:) =  this%i(1,:)

src_mtln/mtln_solver.F90

@@ -70,15 +70,11 @@ subroutine initNodes(this)
         end do
     end subroutine
 
-    subroutine mtln_step(this, currents, voltages)
+    subroutine mtln_step(this)
         class(mtln_t) :: this
-        ! real, dimension(:,:), intent(inout) :: currents
-        ! real, dimension(:,:), intent(in) :: voltages
-        real, dimension(:,:,:), intent(inout) :: currents
-        real, dimension(:,:,:), intent(in) :: voltages
         integer :: i 
 
-        call this%setExternalVoltage(voltages)
+        call this%setExternalVoltage()
 
         call this%advanceBundlesVoltage()
         call this%advanceNWVoltage()
@@ -87,7 +83,7 @@ subroutine mtln_step(this, currents, voltages)
         call this%advanceTime()
         call this%updateProbes()
 
-        call this%updateExternalCurrent(currents)
+        ! call this%updateExternalCurrent()
     end subroutine
 
     subroutine step_alone(this)
@@ -104,14 +100,11 @@ subroutine step_alone(this)
 
     end subroutine
 
-    subroutine setExternalVoltage(this, voltages)
+    subroutine setExternalVoltage(this)
         class(mtln_t) :: this
-        ! real, dimension(:,:), intent(in) :: voltages
-        real, dimension(:,:,:), intent(in) :: voltages
         integer :: i
         do i = 1, this%number_of_bundles
-            call this%bundles(i)%setExternalVoltage(voltages)
-            ! call this%bundles(i)%setExternalVoltage(voltages(i,:))
+            call this%bundles(i)%setExternalVoltage()
         end do
 
     end subroutine

src_mtln/mtln_types.F90

@@ -28,14 +28,14 @@ module mtln_types_mod
    integer, parameter :: DIRECTION_Z_POS   =  3
    integer, parameter :: DIRECTION_Z_NEG   =  -3
 
-   type :: segment_relative_position_t
+   type :: external_field_segment_t
       integer, dimension(3) ::position
       integer :: direction     
-       REAL (KIND=RKIND) , pointer                 ::  Efield_wire2main,Efield_main2wire   !sggmtln
+      real (kind=rkind) , pointer                 ::  Efield_wire2main,Efield_main2wire   !sggmtln
    contains
       private
-      procedure :: segment_relative_positions_eq
-      generic, public :: operator(==) => segment_relative_positions_eq
+      procedure :: external_field_segments_eq 
+      generic, public :: operator(==) => external_field_segments_eq
    end type
 
    type :: termination_t
@@ -120,7 +120,7 @@ module mtln_types_mod
       integer :: conductor_in_parent = -1
       type(connector_t), pointer :: initial_connector => null()
       type(connector_t), pointer :: end_connector => null()
-      type(segment_relative_position_t), allocatable, dimension(:) :: segment_relative_positions
+      type(external_field_segment_t), allocatable, dimension(:) :: external_field_segments
 
    contains
       private
@@ -184,7 +184,7 @@ recursive elemental logical function cable_eq(a,b)
       cable_eq = cable_eq .and.  all(a%step_size == b%step_size)
       cable_eq = cable_eq .and.  (a%transfer_impedance == b%transfer_impedance)
       cable_eq = cable_eq .and.  (a%conductor_in_parent == b%conductor_in_parent)
-      cable_eq = cable_eq .and.  all(a%segment_relative_positions == b%segment_relative_positions)
+      cable_eq = cable_eq .and.  all(a%external_field_segments == b%external_field_segments)
 
 
       if (.not. cable_eq) then
@@ -307,11 +307,31 @@ elemental logical function terminal_network_eq(a,b)
          all(a%connections == b%connections)
    end function
 
-   elemental logical function segment_relative_positions_eq(a,b)
-      class(segment_relative_position_t), intent(in) :: a,b
-      segment_relative_positions_eq = &
+   elemental logical function external_field_segments_eq(a,b)
+      class(external_field_segment_t), intent(in) :: a,b
+      external_field_segments_eq = &
          all(a%position == b%position) .and. &
          a%direction == b%direction
+
+      if (.not. associated(a%Efield_main2wire) .and. .not. associated(b%Efield_main2wire)) then
+         external_field_segments_eq = external_field_segments_eq .and. .true.
+      else if ((associated(a%Efield_main2wire) .and. .not. associated(b%Efield_main2wire)) .or. &
+         (.not. associated(a%Efield_main2wire) .and. associated(b%Efield_main2wire))) then
+            external_field_segments_eq = external_field_segments_eq .and. .false.
+      else
+         external_field_segments_eq = external_field_segments_eq .and. (a%Efield_main2wire == b%Efield_main2wire)
+      end if
+
+      if (.not. associated(a%Efield_wire2main) .and. .not. associated(b%Efield_wire2main)) then
+         external_field_segments_eq = external_field_segments_eq .and. .true.
+      else if ((associated(a%Efield_wire2main) .and. .not. associated(b%Efield_wire2main)) .or. &
+         (.not. associated(a%Efield_wire2main) .and. associated(b%Efield_wire2main))) then
+            external_field_segments_eq = external_field_segments_eq .and. .false.
+      else
+         external_field_segments_eq = external_field_segments_eq .and. (a%Efield_wire2main == b%Efield_wire2main)
+      end if
+
+
    end function
 
    subroutine terminal_connection_add_node(this, node)

src_mtln/preprocess.F90

@@ -211,7 +211,7 @@ function buildLineFromCable(cable) result(res)
                              parent_name = parent_name, &
                              conductor_in_parent = conductor_in_parent, & 
                              transfer_impedance = cable%transfer_impedance, &
-                             segment_relative_positions = cable%segment_relative_positions)
+                             external_field_segments = cable%external_field_segments)
 
         if (associated(cable%initial_connector)) call addConnector(res, cable%initial_connector, 0)
         if (associated(cable%end_connector))     call addConnector(res, cable%initial_connector, size(res%rpul,1))