Merge pull request #249 from OpenSEMBA/dev

Merge dev into master

Author: adrianarce-elemwave

.github/workflows/ubuntu.yml

@@ -46,14 +46,15 @@ jobs:
           #   mtln: "Yes"
           #   hdf: "Yes"
           #   double-precision: "No"
-
-          - os: ubuntu-latest
-            compiler: {name: 'nvidia-hpc', version: '24.5'}
-            build-type: "Release"
-            mpi: "No"
-            mtln: "No"
-            hdf: "No"
-            double-precision: "No"
+          
+          # Disable by lack of space on github action
+          # - os: ubuntu-latest
+          #   compiler: {name: 'nvidia-hpc', version: '24.5'}
+          #   build-type: "Release"
+          #   mpi: "No"
+          #   mtln: "No"
+          #   hdf: "No"
+          #   double-precision: "No"
 
           - os: ubuntu-latest # This is the only test with double precision.
             compiler: {name: 'intel', version: '2025.1'}

CMakeLists.txt

@@ -60,7 +60,7 @@ if (CMAKE_SYSTEM_NAME MATCHES "Linux")
 		set(CMAKE_CXX_FLAGS_RELEASE "-Ofast")
 		set(CMAKE_Fortran_FLAGS_RELEASE "-Ofast")
 
-		set(CMAKE_Fortran_FLAGS_DEBUG "-g -O0")
+		set(CMAKE_Fortran_FLAGS_DEBUG "-g -O0 -fno-inline -fcheck=all -fbacktrace")
 
 	elseif(CMAKE_Fortran_COMPILER_ID MATCHES "IntelLLVM")
 		message(STATUS "Using IntelLLVM (ifx) flags")

src_json_parser/smbjson.F90

@@ -1430,6 +1430,8 @@ function readBlockProbe(bp) result(res)
          type(coords), dimension(:), allocatable :: cs
          type(cell_region_t), dimension(:), allocatable :: cRs
 
+         character(len=1) :: direction
+
          cRs = this%mesh%getCellRegions(this%getIntsAt(bp, J_ELEMENTIDS))
          if (size(cRs) /= 1) then
             call WarnErrReport("Bulk current probe must be defined by a single cell region.", .true.)
@@ -1447,6 +1449,15 @@ function readBlockProbe(bp) result(res)
          res%k1  = cs(1)%zi
          res%k2  = cs(1)%ze
          res%nml = abs(cs(1)%Or)
+         if (res%nml == 0) then !DIR_NULL 
+            direction = this%getStrAt(bp, J_DIR)
+            select case(trim(adjustl(direction)))
+            case(J_DIR_X); res%nml = 1 !DIR_X
+            case(J_DIR_Y); res%nml = 2 !DIR_Y
+            case(J_DIR_Z); res%nml = 3 !DIR_Z
+            case default; call WarnErrReport('Null direction detected for bulk probe. Check definition')
+            end select
+         end if
 
          res%outputrequest = trim(adjustl(this%getStrAt(bp, J_NAME)))
          call setDomain(res, this%getDomain(bp, J_PR_DOMAIN))
@@ -2722,6 +2733,7 @@ function buildNode(termination_list, label, index, id) result(res)
          type(polyline_t) :: polyline
          type(aux_node_t) :: res
          integer :: cable_index
+         integer :: stat
          call this%core%get_child(termination_list, index, termination)
 
          res%node%termination%termination_type = readTerminationType(termination)
@@ -2735,17 +2747,18 @@ function buildNode(termination_list, label, index, id) result(res)
          res%node%side = label
          res%node%conductor_in_cable = index
 
-         call elemIdToCable%get(key(id), value=cable_index)
+         call elemIdToCable%get(key(id), value=cable_index, stat=stat)
+         if (stat == 0) then
          res%node%belongs_to_cable => mtln_res%cables(cable_index)%ptr
 
          polyline = this%mesh%getPolyline(id)
-
-         if (label == TERMINAL_NODE_SIDE_INI) then
-            res%cId = polyline%coordIds(1)
-            res%relPos = this%mesh%getCoordinate(polyline%coordIds(1))
-         else if (label == TERMINAL_NODE_SIDE_END) then
-            res%cId = polyline%coordIds(ubound(polyline%coordIds,1))
-            res%relPos = this%mesh%getCoordinate(polyline%coordIds(ubound(polyline%coordIds,1)))
+            if (label == TERMINAL_NODE_SIDE_INI) then
+               res%cId = polyline%coordIds(1)
+               res%relPos = this%mesh%getCoordinate(polyline%coordIds(1))
+            else if (label == TERMINAL_NODE_SIDE_END) then
+               res%cId = polyline%coordIds(ubound(polyline%coordIds,1))
+               res%relPos = this%mesh%getCoordinate(polyline%coordIds(ubound(polyline%coordIds,1)))
+            end if
          end if
       end function
 
@@ -3176,9 +3189,10 @@ function getPointerToParentCable(cables, id) result(res)
          call elemIdToCable%check_key(key(id), mStat)
          if (mStat /= 0) then
             res => null()
+         else
+            call elemIdToCable%get(key(id), value=index)
+            res => cables(index)%ptr
          end if
-         call elemIdToCable%get(key(id), value=index)
-         res => cables(index)%ptr
       end function
 
       function findConnectorWithId(conn_Id) result(res)
@@ -3480,29 +3494,54 @@ function buildSegments(j_cable, despl) result(res)
          end do
       end function
 
+      pure integer function clip(i, lo, hi)
+         integer, intent(in) :: i, lo, hi
+         clip = max(lo, min(i, hi))
+      end function clip
+
       function getdualBoxYZ(segment, despl) result (res)
          type(Desplazamiento), intent(in) :: despl
          type(segment_t), intent(in) :: segment
          type(box_2d_t) :: res
-         res%min = [-0.5*despl%desY(segment%y-1),-0.5*despl%desZ(segment%z-1)]
-         res%max = [ 0.5*despl%desY(segment%y),   0.5*despl%desZ(segment%z)]
+         integer :: y0, y1, z0, z1
+
+         y0 = clip(segment%y-1, 0, size(despl%desY)-1)
+         y1 = clip(segment%y,   0, size(despl%desY)-1)
+         z0 = clip(segment%z-1, 0, size(despl%desZ)-1)
+         z1 = clip(segment%z,   0, size(despl%desZ)-1)
+
+         res%min = [-0.5 * despl%desY(y0), -0.5 * despl%desZ(z0)]
+         res%max = [ 0.5 * despl%desY(y1),  0.5 * despl%desZ(z1)]
       end function
 
       function getdualBoxXY(segment, despl) result (res)
          type(Desplazamiento), intent(in) :: despl
          type(segment_t), intent(in) :: segment
          type(box_2d_t) :: res
-         res%min = [-0.5*despl%desX(segment%x-1),-0.5*despl%desY(segment%y-1)]
-         res%max = [ 0.5*despl%desX(segment%x),   0.5*despl%desY(segment%y)]
+         integer :: x0, x1, y0, y1
+
+         x0 = clip(segment%x-1, 0, size(despl%desX)-1)
+         x1 = clip(segment%x,   0, size(despl%desX)-1)
+         y0 = clip(segment%y-1, 0, size(despl%desY)-1)
+         y1 = clip(segment%y,   0, size(despl%desY)-1)
+
+         res%min = [-0.5 * despl%desX(x0), -0.5 * despl%desY(y0)]
+         res%max = [ 0.5 * despl%desX(x1),  0.5 * despl%desY(y1)]
       end function
 
       function getdualBoxZX(segment, despl) result (res)
          type(Desplazamiento), intent(in) :: despl
          type(segment_t), intent(in) :: segment
          type(box_2d_t) :: res
-         res%min = [-0.5*despl%desZ(segment%z-1),-0.5*despl%desX(segment%x-1)]
-         res%max = [ 0.5*despl%desZ(segment%z),   0.5*despl%desX(segment%x)]
+         integer :: z0, z1, x0, x1
+
+         z0 = clip(segment%z-1, 0, size(despl%desZ)-1)
+         z1 = clip(segment%z,   0, size(despl%desZ)-1)
+         x0 = clip(segment%x-1, 0, size(despl%desX)-1)
+         x1 = clip(segment%x,   0, size(despl%desX)-1)
 
+         res%min = [-0.5 * despl%desZ(z0), -0.5 * despl%desX(x0)]
+         res%max = [ 0.5 * despl%desZ(z1),  0.5 * despl%desX(x1)]
       end function
 
       function buildStepSize(segments, despl) result(res)

src_json_parser/smbjson_labels.F90

@@ -8,6 +8,7 @@ module smbjson_labels_mod
    character (len=*), parameter :: J_TYPE = "type"
    character (len=*), parameter :: J_ELEMENTIDS = "elementIds"
 
+   character (len=*), parameter :: J_DIR = "direction"
    character (len=*), parameter :: J_DIR_X = "x"
    character (len=*), parameter :: J_DIR_Y = "y"
    character (len=*), parameter :: J_DIR_Z = "z"

test/observation/observation_testingTools.F90

@@ -59,7 +59,11 @@ subroutine check_shape_real(arr, n_expected, test_err, name)
       integer, allocatable :: shp(:)
       character(len=:), allocatable :: nm
 
-      nm = merge(name, "array", present(name))
+      if (present(name)) then
+          nm = trim(adjustl(name))
+      else
+          nm = "array"
+      end if
 
       rank_arr = rank(arr)
       shp = shape(arr)
@@ -82,7 +86,11 @@ subroutine check_shape_complex(arr, n_expected, test_err, name)
       integer, allocatable :: shp(:)
       character(len=:), allocatable :: nm
 
-      nm = merge(name, "array", present(name))
+      if (present(name)) then
+          nm = trim(adjustl(name))
+      else
+          nm = "array"
+      end if
 
       rank_arr = rank(arr)
       shp = shape(arr)
@@ -105,7 +113,11 @@ subroutine check_size(arr, n_expected, test_err, name)
       integer :: siz
       character(len=:), allocatable :: nm
 
-      nm = merge(name, "array", present(name))
+      if (present(name)) then
+          nm = trim(adjustl(name))
+      else
+          nm = "array"
+      end if
 
       rank_arr = rank(arr)
       siz = size(arr)
@@ -200,7 +212,7 @@ function create_base_sgg() result(sgg)
       type(SGGFDTDINFO) :: sgg
 
       sgg%NumMedia = 3
-      allocate(sgg%Med(1:sgg%NumMedia))
+      allocate(sgg%Med(0:sgg%NumMedia))
       sgg%Med = create_basic_media()
       sgg%NumberRequest = 1
       sgg%dt = 0.1_RKIND_tiempo

test/observation/test_observation_init.F90

@@ -29,7 +29,7 @@ integer function test_init_time_movie_observation() bind(C) result(err)
     ThereAreWires = .false.
     ThereAreFarFields = .false.
 
-    initialtimestep = 0
+    initialtimestep = 1
     lastexecutedtime = 0.0_RKIND_tiempo
 
     SINPML_fullsize = create_limit_t(0,4,0,4,0,4,3,3,3)

test/observation/test_observation_update.F90

@@ -31,7 +31,7 @@ integer function test_update_time_movie_observation() bind(C) result(err)
     ThereAreWires = .false.
     ThereAreFarFields = .false.
 
-    initialtimestep = 0
+    initialtimestep = 1
     lastexecutedtime = 0.0_RKIND_tiempo
 
     SINPML_fullsize = create_limit_t(0,4,0,4,0,4,3,3,3)

test/pyWrapper/test_full_system.py

@@ -1024,4 +1024,35 @@ def test_negative_offset_in_x(tmp_path):
 
     assert np.corrcoef(probeTotal['current'].to_numpy(), I_interp)[0, 1] > 0.999
     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)
+    assert np.allclose(probeR['current'].to_numpy(), 0.0, atol=3e-3)
+
+def test_bulk_current_outputs(tmp_path):
+    # This test uses bulk_probe_cases_over_nodal_source.fdtd from input_examples as input.
+    # Verifies all kind of bulk probes are recognised and setted properly by checking outputFile format.
+    fn = PROBES_INPUT_EXAMPLE + 'bulk_probe_cases_over_nodal_source.fdtd.json'
+    solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path)
+    solver.run()
+
+    bulkXPlaneFiles = solver.getSolvedProbeFilenames("BulkXPlane") 
+    bulkYPlaneFiles = solver.getSolvedProbeFilenames("BulkYPlane") 
+    bulkZPlaneFiles = solver.getSolvedProbeFilenames("BulkZPlane") 
+    bulkYPointFiles = solver.getSolvedProbeFilenames("BulkYPoint") 
+    bulkZVolumeFiles = solver.getSolvedProbeFilenames("BulkZVolume") 
+
+    assert len(bulkXPlaneFiles) == 1
+    assert len(bulkYPlaneFiles) == 1
+    assert len(bulkZPlaneFiles) == 1
+    assert len(bulkYPointFiles) == 1
+    assert len(bulkZVolumeFiles) == 10
+
+    probeBulkXPlane = Probe(bulkXPlaneFiles[0])
+    probeBulkYPlane = Probe(bulkYPlaneFiles[0])
+    probeBulkZPlane = Probe(bulkZPlaneFiles[0])
+    probeBulkYPoint = Probe(bulkYPointFiles[0])
+    probeBulkZVolume = Probe(bulkZVolumeFiles[0])
+
+    assert probeBulkXPlane.direction == 'x'
+    assert probeBulkYPlane.direction == 'y'
+    assert probeBulkZPlane.direction == 'z'
+    assert probeBulkYPoint.direction == 'y'
+    assert probeBulkZVolume.direction == 'z'

test/pyWrapper/utils.py

@@ -41,6 +41,7 @@
 OUTPUTS_FOLDER = os.path.join(TEST_DATA_FOLDER, 'outputs/')
 SPINIT_FOLDER = os.path.join(TEST_DATA_FOLDER, 'spinit/')
 GEOMETRIES_FOLDER = os.path.join(TEST_DATA_FOLDER, 'geometries/')
+PROBES_INPUT_EXAMPLE = os.path.join(TEST_DATA_FOLDER, 'input_examples/probes/')
 
 def getCase(case):
     return json.load(open(CASES_FOLDER + case + '.fdtd.json'))

test/smbjson/test_read_mtln.F90

@@ -129,7 +129,7 @@ function expectedProblemDescription() result (expected)
       ! cable 1 - wire
       allocate(unshielded_multiwire_t :: expected%mtln%cables(1)%ptr)
       ptr => expected%mtln%cables(1)%ptr
-      call initializeCablePULParameters(ptr)
+      call initializeCablePULParameters(ptr, n=1)
       select type(ptr)
       type is (unshielded_multiwire_t)
          ptr%name = "line_0_0"
@@ -156,7 +156,7 @@ function expectedProblemDescription() result (expected)
       ! cable 2 - shieldedMultiwire
       allocate(shielded_multiwire_t :: expected%mtln%cables(2)%ptr)
       ptr => expected%mtln%cables(2)%ptr
-      call initializeCablePULParameters(ptr)
+      call initializeCablePULParameters(ptr, n=1)
       select type(ptr)
       type is (shielded_multiwire_t)
          ptr%name = "line_1_0"
@@ -244,7 +244,7 @@ function expectedProblemDescription() result (expected)
       ! cable 4 - wire
       allocate(unshielded_multiwire_t :: expected%mtln%cables(4)%ptr)
       ptr => expected%mtln%cables(4)%ptr
-      call initializeCablePULParameters(ptr)
+      call initializeCablePULParameters(ptr, n=1)
       select type(ptr)
       type is (unshielded_multiwire_t)
          ptr%name = "line_0_1"
@@ -272,7 +272,7 @@ function expectedProblemDescription() result (expected)
       ! cable 5 - shieldedMultiwire
       allocate(shielded_multiwire_t :: expected%mtln%cables(5)%ptr)
       ptr => expected%mtln%cables(5)%ptr
-      call initializeCablePULParameters(ptr)
+      call initializeCablePULParameters(ptr, n=1)
       select type(ptr)
       type is (shielded_multiwire_t)
          ptr%name = "line_1_1"
@@ -346,7 +346,7 @@ function expectedProblemDescription() result (expected)
       ! cable 7 - wire
       allocate(unshielded_multiwire_t :: expected%mtln%cables(7)%ptr)
       ptr => expected%mtln%cables(7)%ptr
-      call initializeCablePULParameters(ptr)
+      call initializeCablePULParameters(ptr, n=1)
       select type(ptr)
       type is (unshielded_multiwire_t)
          ptr%name = "line_0_2"
@@ -374,7 +374,7 @@ function expectedProblemDescription() result (expected)
       ! cable 8 - shieldedMultiwire
       allocate(shielded_multiwire_t :: expected%mtln%cables(8)%ptr)
       ptr => expected%mtln%cables(8)%ptr
-      call initializeCablePULParameters(ptr)
+      call initializeCablePULParameters(ptr, n=1)
       select type(ptr)
       type is (shielded_multiwire_t)
          ptr%name = "line_1_2"
@@ -408,7 +408,7 @@ function expectedProblemDescription() result (expected)
       ! cable 9 - shieldedMultiwire
       allocate(shielded_multiwire_t :: expected%mtln%cables(9)%ptr)
       ptr => expected%mtln%cables(9)%ptr
-      call initializeCablePULParameters(ptr,6)
+      call initializeCablePULParameters(ptr,n=6)
       select type(ptr)
       type is (shielded_multiwire_t)
          ptr%name = "line_2_5"