Implements conductivity for 1D DGTD

.gitignore

@@ -3,4 +3,4 @@
 maxwell/*/__pycache__/*
 test/*/__pycache__/*
 integrators/__pycache__/*
-fdtd/__pycache__
+fdtd/__pycache__/*

maxwell/dg/dg1d.py

@@ -6,7 +6,7 @@
 
 
 class DG1D(SpatialDiscretization):
-    def __init__(self, n_order: int, mesh: Mesh1D, fluxType="Upwind"):
+    def __init__(self, n_order: int, mesh: Mesh1D, fluxType="Upwind",epsilon=None,sigma=None):
         SpatialDiscretization.__init__(self, mesh)
 
         assert n_order > 0
@@ -21,8 +21,24 @@ def __init__(self, n_order: int, mesh: Mesh1D, fluxType="Upwind"):
         alpha = 0
         beta = 0
 
-        # Set up material parameters
-        self.epsilon = np.ones(mesh.number_of_elements())
+        # Epsilon implementation in 1D
+        if epsilon is None:
+            self.epsilon = np.ones(mesh.number_of_elements())
+        elif len(epsilon) != mesh.number_of_elements():
+            raise ValueError("The dimensions of the permittivity vector must align with the number of elements in the mesh.")
+        else:          
+            self.epsilon = np.array(epsilon)
+
+
+        # Sigma implementation necessary for J for 1D
+        if sigma is None:
+            self.sigma = np.zeros(mesh.number_of_elements())
+        elif len(sigma) != mesh.number_of_elements():
+            raise ValueError("The dimensions of the charge density vector must align with the number of elements in the mesh.")
+        else:          
+            self.sigma = np.array(sigma)       
+
+
         self.mu = np.ones(mesh.number_of_elements())
 
         self.x = nodes_coordinates(n_order, mesh.EToV, mesh.vx)
@@ -72,6 +88,7 @@ def buildFields(self):
                       self.mesh.number_of_elements()])
         H = np.zeros(E.shape)
 
+
         return {"E": E, "H": H}
 
     def get_impedance(self):
@@ -152,12 +169,15 @@ def computeJumps(self, E, H):
     def computeRHSE(self, fields):
         E = fields['E']
         H = fields['H']
+        J = np.zeros((self.number_of_nodes_per_element(), self.mesh.number_of_elements()))
+
+        J[:, :] = E * self.sigma
 
         flux_E = self.computeFluxE(E, H)
         rhs_drH = np.matmul(self.diff_matrix, H)
         rhsE = 1/self.epsilon * \
             (np.multiply(-1*self.rx, rhs_drH) +
-             np.matmul(self.lift, self.f_scale * flux_E))
+             np.matmul(self.lift, self.f_scale * flux_E)-J) 
 
         return rhsE
 
@@ -169,8 +189,10 @@ def computeRHSH(self, fields):
         rhs_drE = np.matmul(self.diff_matrix, E)
         rhsH = 1/self.mu * (np.multiply(-1*self.rx, rhs_drE) +
                             np.matmul(self.lift, self.f_scale * flux_H))
+
         return rhsH
 
+
     def computeRHS(self, fields):
         rhsE = self.computeRHSE(fields)
         rhsH = self.computeRHSH(fields)

maxwell/driver.py

@@ -38,7 +38,7 @@ def __init__(self,
 
         # Init time integrator
         if timeIntegratorType == 'EULER':
-            self.timeIntegrator = EULER(self.sp, self.fields)
+            self.timeIntegrator = EULER(self.sp, self.fields)   
         elif timeIntegratorType == 'LSERK4':
             self.timeIntegrator = LSERK4(self.sp, self.fields)
         elif timeIntegratorType == 'LSERK74':

requirements.txt

@@ -2,4 +2,5 @@ pytest >= 7.3
 numpy
 matplotlib
 scipy
-nodepy
+nodepy
+scikit-rf

test/test_DFT.py

@@ -0,0 +1,64 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+import pytest
+import time
+
+def dft(x):
+    N = len(x)
+    X = np.zeros(N, dtype=complex)
+    for i in range(N):
+        for n in range(N):
+            X[i]=X[i]+ x[n] * np.exp(-2j*np.pi*i*n/N)
+    return X
+
+def test_DFT():
+    # Function to manually compute the DFT
+
+    # Create the Gaussian function
+    x = np.linspace(-5, 5, 100) 
+    s0 = 0.50
+    gaussian = np.exp(-(x) ** 2 / (2 * s0 **2))  
+
+    # Compue FFT and DFT
+    fft_g = np.fft.fft(gaussian)
+
+    dft_g = dft(gaussian)
+
+    # Assert to check the difference between both transforms
+    assert np.allclose(dft_g, fft_g, atol=0.01), "DFT and FFT differ by more than 0.01"
+
+    # # Graph
+    # plt.figure(figsize=(12, 6))
+    # # Ejes de frecuencia para FFT y DFT
+    # freq = np.fft.fftshift(np.fft.fftfreq(len(x), x[1] - x[0]))
+
+    # # Original
+    # plt.subplot(1, 3, 1)
+    # plt.plot(x, gaussian, label='Gaussian')
+    # plt.title('Gaussian')
+    # plt.xlabel('x')
+    # plt.ylabel('Amplitude')
+    # plt.grid()
+    # plt.legend()
+
+    # # DFT
+    # plt.subplot(1, 3, 2)
+    # plt.plot(freq, np.abs(fft_g), label='FFT')
+    # plt.title('FFT Transform')
+    # plt.xlabel('Frequency')
+    # plt.ylabel('Amplitude')
+    # plt.grid()
+    # plt.legend()
+
+    # # FFT
+    # plt.subplot(1, 3, 3)
+    # plt.plot(freq, np.abs(dft_g), label='DFT', linestyle='--')
+    # plt.title('DFT Transform')
+    # plt.xlabel('Frequency')
+    # plt.ylabel('Amplitude')
+    # plt.grid()
+    # plt.legend()
+
+    # plt.tight_layout()
+    # plt.show()