266 Wires tangential to surfaces

CMakePresets.json

@@ -24,6 +24,14 @@
       "cacheVariables": {
         "SEMBA_FDTD_ENABLE_MTLN": "OFF"
       }
+    },
+    {
+      "name": "rls-nomtln",
+      "inherits": "rls",
+      "binaryDir": "build-rls-nomtln/",
+      "cacheVariables": {
+        "SEMBA_FDTD_ENABLE_MTLN": "OFF"
+      }
     }
   ]
 }

src_main_pub/timestepping.F90

@@ -1,19 +1,3 @@
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!  SEMBA_FDTD sOLVER module
-!  Creation date Date :  April, 8, 2010
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!__________________________________________________________________________________________________
-!******************************** REVISAR PARA PGI (CRAY) *****************************************
-!---> AdvanceMultiportE
-!---> AdvanceAnisMultiportE
-!---> AdvanceMultiportH
-!---> AdvanceAnisMultiportH
-!---> dfUpdateE
-!---> dfUpdateH
-!---> MinusCloneMagneticPMC
-!________________________________________________________________________________________
-
 module Solver_mod
 
    use fdetypes

src_pyWrapper/pyWrapper.py

@@ -119,11 +119,11 @@ def __init__(self, probe_filename):
             self.cell_init, self.cell_end = \
                 Probe._positionStrToTwoCells(position_str)
             self.data = self.data.rename(columns={
-                    self.data.columns[0]: 'freq',
-                    self.data.columns[7]: 'rcs_arit',
-                    self.data.columns[8]: 'rcs_geom'
-                })
-                
+                self.data.columns[0]: 'freq',
+                self.data.columns[7]: 'rcs_arit',
+                self.data.columns[8]: 'rcs_geom'
+            })
+
         elif tag in Probe.MOVIE_TAGS:
             self.type = 'movie'
             self.cell_init, self.cell_end = \
@@ -147,6 +147,11 @@ def __init__(self, probe_filename):
             raise ValueError("Unable to determine probe type")
 
     def __getitem__(self, key):
+        if key not in self.data.columns:
+            if key == 'current_0' and 'current' in self.data.columns:
+                return self.data['current']
+            if key == 'voltage_0' and 'voltage' in self.data.columns:
+                return self.data['voltage']
         return self.data[key]
 
     @staticmethod
@@ -207,6 +212,7 @@ def _positionStrToTwoCells(pos_str: str):
     def _getFieldAndDirection(tag: str):
         return tag[1], tag[2]
 
+
 class ExcitationFile():
     def __init__(self, excitation_filename):
         if isinstance(excitation_filename, os.PathLike):
@@ -215,11 +221,12 @@ def __init__(self, excitation_filename):
             self.filename = excitation_filename
         assert os.path.isfile(self.filename)
 
-        self.data = pd.read_csv(self.filename, sep='\\s+', names=['time', 'value'])
+        self.data = pd.read_csv(
+            self.filename, sep='\\s+', names=['time', 'value'])
 
     def __getitem__(self, key):
         return self.data[key]
-    
+
 
 class FDTD():
     def __init__(self, input_filename, path_to_exe=None,
@@ -351,7 +358,6 @@ def cleanUp(self):
         for f in subfolders:
             shutil.rmtree(f, ignore_errors=True)
 
-
     def getSolvedProbeFilenames(self, probe_name):
         if not "probes" in self._input:
             raise ValueError('Solver does not contain probes.')
@@ -366,18 +372,19 @@ def getSolvedProbeFilenames(self, probe_name):
         return sorted(probeFiles)
 
     def getExcitationFile(self, excitation_file_name):
-        file_extensions =('*.1.exc',)
+        file_extensions = ('*.1.exc',)
         excitationFile = []
         for ext in file_extensions:
-            newExcitationFile = [x for x in glob.glob(ext) if re.match(excitation_file_name, x)]
+            newExcitationFile = [x for x in glob.glob(
+                ext) if re.match(excitation_file_name, x)]
             excitationFile.extend(newExcitationFile)
-        
+
         if ((len(excitationFile)) != 1):
-            raise "Unexpected number of excitation Files found: {}".format(excitationFile)
+            raise ValueError(
+                "Unexpected number of excitation Files found: {}".format(excitationFile))
 
         return excitationFile
 
-
     def getVTKMap(self):
         current_path = os.getcwd()
         folders = [item for item in os.listdir(
@@ -387,18 +394,18 @@ def getVTKMap(self):
         for folder in folders:
             mapFile = os.path.join(current_path, folder, folder+"_1.vtk")
             if os.path.isfile(mapFile):
-                return mapFile            
+                return mapFile
         raise ValueError("Unable to find mapvatk file")
 
-
     def getCurrentVTKMap(self):
         current_path = os.getcwd()
         folders = [item for item in os.listdir(
             current_path) if os.path.isdir(os.path.join(current_path, item))]
         if len(folders) != 1:
             return None
         for folder in folders:
-            mapFile = os.path.join(current_path, folder, folder+"_1_current.vtk")
+            mapFile = os.path.join(current_path, folder,
+                                   folder+"_1_current.vtk")
             if os.path.isfile(mapFile):
                 return mapFile
         raise ValueError("Unable to find current vtk file")
@@ -408,5 +415,3 @@ def getMaterialProperties(self, materialName):
             for idx, element in enumerate(self._input['materials']):
                 if element.get("name") == materialName:
                     return self._input['materials'][idx]
-
-

test/pyWrapper/test_full_system.py

@@ -405,6 +405,68 @@ def test_towelHanger(tmp_path):
     assert np.allclose(p_expected[2]['current'],  p_solved[2]['current_0'], rtol=5e-3, atol=5e-4)
 
 
+def test_towel_rack_with_and_without_shorting_plane(tmp_path):
+    fn = CASES_FOLDER + \
+        'towel_rack_with_shorting_plane/towel_rack_with_shorting_plane.fdtd.json'
+
+    # --- excitation ---
+    dt = 1e-12
+    w0 = 0.1e-8
+    t0 = 10 * w0
+    t = np.arange(0, t0 + 20 * w0, dt)
+    excitation = np.exp(-np.power(t - t0, 2) / w0 ** 2)
+    exc_file = os.path.join(tmp_path, 'gauss.exc')
+    np.savetxt(exc_file, np.column_stack([t, excitation]))
+
+    freqs = np.geomspace(1e3, 1e9, 61)
+
+    # --- with shorting plane ---
+    folder_with    = os.path.join(tmp_path, 'with_shorting_plane')
+    os.makedirs(folder_with)
+    solver_w = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=folder_with)
+    solver_w.run()
+
+    probe_name = os.path.basename(solver_w.getSolvedProbeFilenames("Wire probe")[0])
+    probe_w = Probe(os.path.join(folder_with, probe_name))
+    time_I_w  = probe_w["time"].to_numpy()
+    current_w = probe_w["current_0"].to_numpy()
+    V_interp_w = np.interp(time_I_w, t, excitation)
+    Z_in_w = dtft(V_interp_w, time_I_w, freqs) / dtft(current_w, time_I_w, freqs)
+
+    # --- without shorting plane ---
+    folder_without = os.path.join(tmp_path, 'without_shorting_plane')
+    os.makedirs(folder_without)
+    solver_wo = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE, run_in_folder=folder_without)
+    solver_wo['materialAssociations'][0]['elementIds'] = [1]
+    solver_wo.run()
+
+    probe_wo = Probe(os.path.join(folder_without, probe_name))
+    time_I_wo  = probe_wo["time"].to_numpy()
+    current_wo = probe_wo["current_0"].to_numpy()
+    V_interp_wo = np.interp(time_I_wo, t, excitation)
+    Z_in_wo = dtft(V_interp_wo, time_I_wo, freqs) / dtft(current_wo, time_I_wo, freqs)
+
+    # For debugging
+    # import matplotlib.pyplot as plt
+    # plt.figure()
+    # plt.semilogx(freqs, 20 * np.log10(np.abs(Z_in_w)),  label="With shorting plane")
+    # plt.semilogx(freqs, 20 * np.log10(np.abs(Z_in_wo)), label="Without shorting plane")
+    # plt.xlabel("Frequency [Hz]")
+    # plt.ylabel("|Z(j2πf)| [dB]")
+    # plt.grid(which="both")
+    # plt.legend()
+    # plt.tight_layout()
+    # plt.savefig('tmp-plot.png')
+    # plt.show()
+
+    # Expect the shorting plane not chaning the impedance at low frequencies.
+    assert np.allclose(
+        np.abs(Z_in_w[freqs < 1e6]), 
+        np.abs(Z_in_wo[freqs < 1e6]), 
+        rtol=0.1
+    )
+
+
 @no_hdf_skip
 @pytest.mark.hdf
 def test_sphere(tmp_path):

test/pyWrapper/utils.py

@@ -38,7 +38,8 @@
 elif platform == "win32":
     SEMBA_EXE = os.path.join(os.getcwd(), 'build', 'bin', 'semba-fdtd.exe')
 
-NGSPICE_DLL = os.path.join(os.getcwd(), 'precompiled_libraries', 'windows-intel', 'ngspice', 'ngspice.dll')
+NGSPICE_DLL = os.path.join(
+    os.getcwd(), 'precompiled_libraries', 'windows-intel', 'ngspice', 'ngspice.dll')
 TEST_DATA_FOLDER = os.path.join(os.getcwd(), 'testData/')
 CASES_FOLDER = os.path.join(TEST_DATA_FOLDER, 'cases/')
 MODELS_FOLDER = os.path.join(TEST_DATA_FOLDER, 'models/')
@@ -48,6 +49,7 @@
 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'))
 
@@ -103,31 +105,33 @@ def readSpiceFile(spice_file):
             val = np.append(val, float(l.split()[1]))
     return t, val
 
-def createWire(id, r, rpul = 0.0, lpul=0.0):
-    return {"id":id,
+
+def createWire(id, r, rpul=0.0, lpul=0.0):
+    return {"id": id,
             "type": "wire",
-            "radius": r, 
-            "resistancePerMeter": rpul, 
+            "radius": r,
+            "resistancePerMeter": rpul,
             "inductancePerMeter": lpul
             }
 
-def createUnshieldedWire(id, lpul, cpul, rpul = 0.0, gpul = 0.0):
-    return {         
+
+def createUnshieldedWire(id, lpul, cpul, rpul=0.0, gpul=0.0):
+    return {
         "id": id,
         "type": "unshieldedMultiwire",
-        "inductancePerMeter" :  [[lpul]],
-        "capacitancePerMeter" : [[cpul]],
+        "inductancePerMeter":  [[lpul]],
+        "capacitancePerMeter": [[cpul]],
         "resistancePerMeter": [rpul],
         "conductancePerMeter": [gpul]
-        }
+    }
 
 
-def createPropertyDictionary(vtkfile, celltype:int, property:str):
+def createPropertyDictionary(vtkfile, celltype: int, property: str):
     ugrid = pv.UnstructuredGrid(vtkfile)
     objs = np.argwhere(ugrid.celltypes == celltype)  # [i][0]
     if len(objs) == 0:
         return dict()
-    
+
     props = np.array([])
     prop_array = ugrid.cell_data[property]
     for i in range(objs[0][0], objs[-1][0]+1):
@@ -159,21 +163,43 @@ def copyXSpiceModels(temp_dir, sys_name):
     makeCopy(temp_dir, SPINIT_FOLDER + sys_name + 'xtradev.cm')
     makeCopy(temp_dir, SPINIT_FOLDER + sys_name + 'xtraevt.cm')
 
+
 def RCS(fspace, radius):
 
-    def h_sph(n,z):
-        return np.sqrt(0.5*np.pi*z)*h(n+0.5,z)
+    def h_sph(n, z):
+        return np.sqrt(0.5*np.pi*z)*h(n+0.5, z)
 
-    def h_sph_der(n,z):
-        return 0.25*np.pi*(0.5*np.pi*z)**(-0.5)*h(n+0.5,z) + np.sqrt(0.5*np.pi*z)*hp(n+0.5,z)
+    def h_sph_der(n, z):
+        return 0.25*np.pi*(0.5*np.pi*z)**(-0.5)*h(n+0.5, z) + np.sqrt(0.5*np.pi*z)*hp(n+0.5, z)
 
-    a = radius #radius
+    a = radius  # radius
     f = fspace
     l = 3.0e8/f
     beta = 2*np.pi*f/3.0e8
 
     res = 0.0
-    for n in range(1,50):
-        res = res + (-1)**n*(2*n+1)/(h_sph_der(n,beta*a)*h_sph(n,beta*a))
+    for n in range(1, 50):
+        res = res + (-1)**n*(2*n+1)/(h_sph_der(n, beta*a)*h_sph(n, beta*a))
     res = np.abs(res)**2*(l**2/(4*np.pi))
-    return res
+    return res
+
+
+def dtft(signal, time, freqs):
+    """Continuous-time DTFT approximation using trapezoidal integration."""
+    signal = np.asarray(signal)
+    time = np.asarray(time)
+    res = np.empty_like(freqs, dtype=complex)
+    for i, f in enumerate(freqs):
+        res[i] = np.trapz(signal * np.exp(-1j * 2 * np.pi * f * time), time)
+    return res
+
+
+def compute_impedance(probe_path, time_exc, voltage_exc, freqs):
+    probe = Probe(probe_path)
+    time_I = probe["time"].to_numpy()
+    current = probe["current_0"].to_numpy()
+    V_interp = np.interp(time_I, time_exc, voltage_exc)
+    I_f = dtft(current, time_I, freqs)
+    V_f = dtft(V_interp, time_I, freqs)
+    Z = V_f / I_f
+    return time_I, current, Z