diff --git a/.gitignore b/.gitignore index b168ec8..398fc63 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,4 @@ maxwell/*/__pycache__/* test/*/__pycache__/* integrators/__pycache__/* -fdtd/__pycache__ +fdtd/__pycache__/* \ No newline at end of file diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index 44cb17c..78ee226 100644 --- a/maxwell/dg/dg1d.py +++ b/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) diff --git a/maxwell/driver.py b/maxwell/driver.py index 3c9cd0d..f701d2e 100644 --- a/maxwell/driver.py +++ b/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': diff --git a/requirements.txt b/requirements.txt index 42012b6..91b9630 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,5 @@ pytest >= 7.3 numpy matplotlib scipy -nodepy \ No newline at end of file +nodepy +scikit-rf \ No newline at end of file diff --git a/test/test_DFT.py b/test/test_DFT.py new file mode 100644 index 0000000..5d2e226 --- /dev/null +++ b/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() diff --git a/test/test_TFG_1single_slab.py b/test/test_TFG_1single_slab.py new file mode 100644 index 0000000..d6c8ea3 --- /dev/null +++ b/test/test_TFG_1single_slab.py @@ -0,0 +1,470 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import pytest +import time + + +from maxwell.driver import * +from maxwell.dg.mesh1d import * +from maxwell.dg.dg1d import * +from maxwell.fd.fd1d import * + + +from nodepy import runge_kutta_method as rk + + +from skrf import Frequency +from skrf.media import Freespace + + + +def dft(x,freq,time): + X=[] + for f in range(len(freq)): + summatory=0. + for t in range(len(time)): + summatory=summatory + x[t] * np.exp(-2j*np.pi*freq[f]*time[t]) + X.append(summatory) + return X + + + +def test_TFG_ep50_rho1_slab_1cm(): + + # Material distribution + epsilon_1 = 1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + epsilon_r_material = 50.0 + rho_material = 1.0 + + # Mesh + L1 = 0.0 + L2 = 1.0 + elements = 100 + epsilons = epsilon_1*np.ones(elements) + epsilons[49] = epsilon_r_material + sigmas = np.zeros(elements) + sigmas[49] = Z_0/rho_material + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.25 + final_time = 4.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][60]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.thru() ** mat.line(1, unit = 'cm') ** air.thru() + R_ref=20*np.log10(abs(slab.s[:, 0, 0])) + T_ref=20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 1, \varepsilon = 50$, slab width = 1 cm") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + + + +def test_TFG_ep20_rho5_slab_1cm(): + + # Material distribution + epsilon_1 = 1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + epsilon_r_material = 20.0 + rho_material = 5.0 + + # Mesh + L1 = 0.0 + L2 = 1.0 + elements = 100 + epsilons = epsilon_1*np.ones(elements) + epsilons[49] = epsilon_r_material + sigmas = np.zeros(elements) + sigmas[49] = Z_0/rho_material + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.25 + final_time = 4.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][60]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.thru() ** mat.line(1, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 5, \varepsilon = 20$, slab width = 1 cm") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + + + +def test_TFG_ep20_rho5_slab_6cm(): + + # Material distribution + epsilon_1 = 1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + epsilon_r_material = 20.0 + rho_material = 5.0 + + # Mesh + L1=0.0 + L2=1.0 + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[47:53]=epsilon_r_material + sigmas = np.zeros(elements) + sigmas[47:53] = Z_0/rho_material + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.25 + final_time = 4.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][60]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.thru() ** mat.line(6, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 5, \varepsilon = 20$, slab width = 6 cm") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + + + +def test_TFG_ep6_rho8_slab_6cm(): + + # Material distribution + epsilon_1 = 1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + epsilon_r_material = 6.0 + rho_material = 8.0 + + # Mesh + L1=0.0 + L2=1.0 + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[47:53]=epsilon_r_material + sigmas = np.zeros(elements) + sigmas[47:53] = Z_0/rho_material + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.25 + final_time = 4.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][60]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.thru() ** mat.line(6, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 8, \varepsilon = 6$, slab width = 6 cm") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() \ No newline at end of file diff --git a/test/test_TFG_3multi_slab.py b/test/test_TFG_3multi_slab.py new file mode 100644 index 0000000..7522f15 --- /dev/null +++ b/test/test_TFG_3multi_slab.py @@ -0,0 +1,505 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import pytest +import time + + +from maxwell.driver import * +from maxwell.dg.mesh1d import * +from maxwell.dg.dg1d import * +from maxwell.fd.fd1d import * + + +from nodepy import runge_kutta_method as rk + +from skrf import Frequency +from skrf.media import Freespace + + + +def dft(x,freq,time): + X=[] + for f in range(len(freq)): + summatory=0. + for t in range(len(time)): + summatory=summatory + x[t] * np.exp(-2j*np.pi*freq[f]*time[t]) + X.append(summatory) + return X + + + +def test_TFG_ep50_rho1_multislab_1cm(): + + # Material distribution + epsilon_1=1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + + epsilon_r_material = 50 + rho_material = 1 + + # Mesh + L1 = 0.0 + L2 = 1.0 + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[25] = epsilon_r_material + epsilons[50] = epsilons[25] + epsilons[75] = epsilons[25] + sigmas = np.zeros(elements) + sigmas[25] = Z_0/rho_material + sigmas[50] = sigmas[25] + sigmas[75] = sigmas[25] + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.15 + final_time = 25.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + additional_freq_1 = np.linspace(612e6,730e6,40) + freq_vector = np.union1d(freq_vector, additional_freq_1) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][95]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and Reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + + dB_T = 20*np.log10(T) + dB_R = 20*np.log10(R) + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.line(24, unit = 'cm') ** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 1, \varepsilon = 50$, slab width = 1 cm, N=3") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + + + +def test_TFG_ep112_rho7_multislab_1cm(): + + # Material distribution + epsilon_1=1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + + epsilon_r_material = 11.2 + rho_material = 7.0 + + # Mesh + L1 = 0.0 + L2 = 1.0 + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[25] = epsilon_r_material + epsilons[50] = epsilons[25] + epsilons[75] = epsilons[25] + sigmas = np.zeros(elements) + sigmas[25] = Z_0/rho_material + sigmas[50] = sigmas[25] + sigmas[75] = sigmas[25] + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.15 + final_time = 25.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + additional_freq_1 = np.linspace(686e6,706e6,25) + additional_freq_2 = np.linspace(817e6,843e6,25) + freq_vector = np.union1d(freq_vector, additional_freq_1) + freq_vector = np.union1d(freq_vector, additional_freq_2) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][95]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and Reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.line(24, unit = 'cm') ** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 7, \varepsilon = 11.2$, slab width = 1 cm, N=3") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + + + +def test_TFG_ep20_rho5_multislab_1cm(): + + # Material distribution + epsilon_1=1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + + epsilon_r_material = 20.0 + rho_material = 5.0 + + # Mesh + L1 = 0.0 + L2 = 1.0 + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[25] = epsilon_r_material + epsilons[50] = epsilons[25] + epsilons[75] = epsilons[25] + sigmas = np.zeros(elements) + sigmas[25] = Z_0/rho_material + sigmas[50] = sigmas[25] + sigmas[75] = sigmas[25] + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.15 + final_time = 25.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + additional_freq_1 = np.linspace(650e6,675e6,25) + additional_freq_2 = np.linspace(755e6,778e6,25) + freq_vector = np.union1d(freq_vector, additional_freq_1) + freq_vector = np.union1d(freq_vector, additional_freq_2) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + + # DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][95]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + # Transmission and Reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.line(24, unit = 'cm') ** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 5, \varepsilon = 20$, slab width = 1 cm, N=3") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + + + +def test_TFG_ep5_rho15_multislab_1cm(): + + # Material distribution + epsilon_1=1.0 + mu_0 = 4.0*np.pi*1e-7 + eps_0 = 8.854187817e-12 + Z_0 = np.sqrt(mu_0/eps_0) + + + epsilon_r_material = 5.0 + rho_material = 15.0 + + # Mesh + L1 = 0.0 + L2 = 1.0 + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[25] = epsilon_r_material + epsilons[50] = epsilons[25] + epsilons[75] = epsilons[25] + sigmas = np.zeros(elements) + sigmas[25] = Z_0/rho_material + sigmas[50] = sigmas[25] + sigmas[75] = sigmas[25] + + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="SMA"), + epsilon=epsilons, + sigma=sigmas + ) + driver = MaxwellDriver(sp) + + # Type of wave + s0 = 0.025 + x0 = 0.15 + final_time = 25.0 + steps = int(np.ceil(final_time/driver.dt)) + freq_vector = np.logspace(8, 9, 301) + additional_freq_1 = np.linspace(733e6,751e6,25) + additional_freq_2 = np.linspace(908e6,927e6,25) + freq_vector = np.union1d(freq_vector, additional_freq_1) + freq_vector = np.union1d(freq_vector, additional_freq_2) + + + initialFieldE = np.exp(-(sp.x-x0)**2.0/(2.0*s0**2.0)) + initialFieldH = initialFieldE + + # Driver operates + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + E_vector_R = [] + E_vector_T = [] + E_vector_0 = [] + + + # #DFT calculations + time_vector_coeffs = np.linspace(0.0, final_time, steps) + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + for _ in range(steps): + driver.step() + E_vector_R.append(driver['E'][3][5]) + E_vector_T.append(driver['E'][3][95]) + + for t in time_vector_coeffs: + E_vector_0.append(np.exp(-(t-x0)**2.0/(2.0*s0**2.0))) + + time_vector_coeffs_corrected = time_vector_coeffs / 299792458 + + dft_E_R=dft(E_vector_R,freq_vector,time_vector_coeffs_corrected) + dft_E_T=dft(E_vector_T,freq_vector,time_vector_coeffs_corrected) + dft_0=dft(E_vector_0,freq_vector,time_vector_coeffs_corrected) + + #Transmission and Reflection coefficients + T = np.abs(dft_E_T) / np.abs(dft_0) + R = np.abs(dft_E_R) / np.abs(dft_0) + dB_T=20*np.log10(T) + dB_R=20*np.log10(R) + + # Sci-kit simulation + freq_ref = Frequency(100, 1000, 300, 'MHz') + air = Freespace(freq_ref) + mat = Freespace(freq_ref, ep_r = epsilon_r_material, rho = rho_material) + slab = air.line(24, unit = 'cm') ** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.line(24, unit = 'cm')** mat.line(1, unit = 'cm') ** air.thru() + R_ref = 20*np.log10(abs(slab.s[:, 0, 0])) + T_ref = 20*np.log10(abs(slab.s[:, 1, 0])) + freq_ref_values = freq_ref.f_scaled *1e6 + + + # Interpolate computed R and T to match reference frequencies + R_interp = np.interp(freq_ref_values, freq_vector, dB_R) + T_interp = np.interp(freq_ref_values, freq_vector, dB_T) + + + # Assert conditions + assert np.all(np.abs(R_interp - R_ref) <= 0.1), "Computed R deviates too much!" + assert np.all(np.abs(T_interp - T_ref) <= 0.1), "Computed T deviates too much!" + + + # Generate plot + # plt.figure() + # plt.plot(freq_vector / 1e6, dB_T, label='DGTD T', color='blue') + # plt.plot(freq_vector / 1e6, dB_R, label='DGTD R', color='cyan', linestyle='solid') + # plt.plot(freq_ref_values / 1e6, T_ref, label='Sci-kit T', color='magenta', linestyle='dashed') + # plt.plot(freq_ref_values / 1e6, R_ref, label='Sci-kit R', color='red', linestyle='dashed') + # plt.xlabel("Frequency (MHz)") + # plt.ylabel("Magnitude (dB)") + # plt.title(r"$\rho = 15, \varepsilon = 5$, slab width = 1 cm, N=3") + # plt.legend() + # plt.xlim(100, 1000) + # plt.grid() + # plt.show() + diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 4a9555f..e849eb0 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -23,36 +23,179 @@ def test_spatial_discretization_lift(): assert np.allclose(surface_integral_dg(1, jacobiGL(0.0, 0.0, 1)), np.array([[2.0, -1.0], [-1.0, 2.0]])) -def test_pec(): +def test_pec_dielectric_upwind(): + epsilon_1=1 + epsilon_2=2 + mu_1=1 + mu_2=1 + z_1=np.sqrt(mu_1/epsilon_1) + z_2=np.sqrt(mu_2/epsilon_2) + epsilons = epsilon_1*np.ones(100) + epsilons[50:99]=epsilon_2 + sp = DG1D( - n_order=5, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC") + n_order=3, + mesh=Mesh1D(-5.0, 5.0, 100, boundary_label="PEC"), + epsilon=epsilons ) driver = MaxwellDriver(sp) - final_time = 1.999 - s0 = 0.25 - initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) + + #T and R coeffs + T_E=2*z_2/(z_1+z_2) + R_E=(z_2-z_1)/(z_2+z_1) + + final_time = 6 + s0 = 0.50 + initialFieldE = np.exp(-(sp.x+2)**2/(2*s0**2)) + initialFieldH = initialFieldE driver['E'][:] = initialFieldE[:] - finalFieldE = driver['E'] + driver['H'][:] = initialFieldH[:] + + driver.run_until(final_time) + + max_electric_field_2 = np.max(driver['E'][2][50:99]) + assert np.isclose(max_electric_field_2, T_E, atol=0.1) + + min_electric_field_1 = np.min(driver['E'][2][0:49]) + assert np.isclose(min_electric_field_1, R_E, atol=0.1) + # driver['E'][:] = initialFieldE[:] + # driver['H'][:] = initialFieldH[:] + # for _ in range(300): + # driver.step() + # plt.plot(sp.x, driver['E'],'b') + # plt.plot(sp.x, driver['H'],'r') + # plt.ylim(-1, 1.5) + # plt.grid(which='both') + # plt.pause(0.01) + # plt.cla() + +def test_pec_dielectric_upwind_v(): + epsilon_1=1 + epsilon_2=2 + mu_1=1 + mu_2=1 + z_1=np.sqrt(mu_1/epsilon_1) + z_2=np.sqrt(mu_2/epsilon_2) + v_1=1/np.sqrt(epsilon_1*mu_1) + v_2=1/np.sqrt(epsilon_2*mu_2) + L1=-5.0 + L2=5.0 + Lt=abs(L1)+abs(L2) + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[int(elements/2):elements-1]=epsilon_2 + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="PEC"), + epsilon=epsilons, + sigma=np.zeros(100) + ) + driver = MaxwellDriver(sp) + + + #T and R coeffs + T_E=2*z_2/(z_1+z_2) + R_E=(z_2-z_1)/(z_2+z_1) + + s0 = 0.50 + x0=-2 + t_0=abs(x0/v_1) + final_time = 6 + + initialFieldE = np.exp(-(sp.x-x0)**2/(2*s0**2)) + initialFieldH = initialFieldE + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + driver.run_until(final_time) - R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), - -finalFieldE.reshape(1, finalFieldE.size)) - assert R[0, 1] > 0.9999 + reflectedelement=int(abs(v_1*final_time+L1)*elements/Lt) + transelement=int((-L1+(final_time-t_0)*v_2)*elements/Lt) + + max_electric_field_2 = np.max(driver['E'][2][transelement-1:transelement+1]) + assert np.isclose(max_electric_field_2, T_E, atol=0.1) + + min_electric_field_1 = np.min(driver['E'][2][reflectedelement-1:reflectedelement+1]) + assert np.isclose(min_electric_field_1, R_E, atol=0.1) # driver['E'][:] = initialFieldE[:] - # for _ in range(172): + # driver['H'][:] = initialFieldH[:] + # for _ in range(300): # driver.step() # plt.plot(sp.x, driver['E'],'b') # plt.plot(sp.x, driver['H'],'r') - # plt.ylim(-1, 1) + # plt.ylim(-1, 1.5) # plt.grid(which='both') # plt.pause(0.01) # plt.cla() +def test_pec_dielectric_centered_right(): + epsilon_1=2 + epsilon_2=1 + mu_1=1 + mu_2=1 + z_1=np.sqrt(mu_1/epsilon_1) + z_2=np.sqrt(mu_2/epsilon_2) + v_1=1/np.sqrt(epsilon_1*mu_1) + v_2=1/np.sqrt(epsilon_2*mu_2) + L1=-5.0 + L2=5.0 + Lt=abs(L1)+abs(L2) + elements=100 + epsilons = epsilon_1*np.ones(elements) + epsilons[int(elements/2):elements-1]=epsilon_2 + + sp = DG1D( + n_order=3, + mesh=Mesh1D(L1, L2, elements, boundary_label="PEC"), + fluxType="Centered", + epsilon=epsilons + ) + driver = MaxwellDriver(sp) + + + #T and R coeffs + T_E=-2*z_1/(z_1+z_2) + R_E=-(z_1-z_2)/(z_2+z_1) + + s0 = 0.5 + x0=2 + t_front=(2*L2-x0)/v_2 + final_time = 10 + + initialFieldE = np.exp(-(sp.x-x0)**2/(2*s0**2)) + initialFieldH = initialFieldE + + driver['E'][:] = initialFieldE[:] + driver['H'][:] = initialFieldH[:] + + driver.run_until(final_time) + + reflectedelement=int((-L1+v_2*(final_time-t_front))*elements/Lt) + transelement=int((-L1-v_1*(final_time-t_front))*elements/Lt) + + max_electric_field_2 = np.max(driver['E'][2][reflectedelement-1:reflectedelement+1]) + assert np.isclose(max_electric_field_2, R_E, atol=0.01) + + min_electric_field_1 = np.min(driver['E'][2][transelement-1:transelement+1]) + assert np.isclose(min_electric_field_1, T_E, atol=0.01) + + # driver['E'][:] = initialFieldE[:] + # driver['H'][:] = initialFieldH[:] + + # for _ in range(240): + # driver.step() + # plt.plot(sp.x, driver['E'],'b') + # plt.plot(sp.x, driver['H'],'r') + # plt.ylim(-1.5, 1.5) + # plt.grid(which='both') + # plt.pause(0.01) + # plt.cla() def test_pec_centered(): sp = DG1D( @@ -814,7 +957,7 @@ def real_function(x, t): print("ERROR LSERK134", error_RMSE_LSERK134) # Show dt in driver - assert (np.sqrt(error_abs).max() < 1e-02, True) + assert np.max(np.sqrt(error_abs)) < 1e-02 @pytest.mark.skip(reason="Nothing is being tested.")