From d25104bc0c374e964cee879f1847db29bbebe264 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 30 Apr 2024 09:37:11 +0200 Subject: [PATCH 01/25] Adds option to reduce essential DoF in drived operator --- maxwell/driver.py | 11 ++- maxwell/fd/fd1d.py | 154 ++++++++++++++------------------ test/fd/test_fd1d.py | 2 +- test/test_dgtd_1d.py | 2 +- test/tests_fdtd/test_fdtd_1d.py | 19 +++- 5 files changed, 96 insertions(+), 92 deletions(-) diff --git a/maxwell/driver.py b/maxwell/driver.py index 3c9cd0d..26deae4 100644 --- a/maxwell/driver.py +++ b/maxwell/driver.py @@ -12,6 +12,7 @@ from .integrators.LF2V import * from .integrators.EULER import * +import copy class MaxwellDriver: def __init__(self, @@ -79,9 +80,12 @@ def run_until(self, final_time): def __getitem__(self, key): return self.fields[key] - def buildDrivedEvolutionOperator(self): + def buildDrivedEvolutionOperator(self, reduceToEsentialDoF=True): N = self.sp.number_of_unknowns() A = np.zeros((N,N)) + + oldFields = copy.deepcopy(self.fields) + for i in range(N): self.fields = self.sp.buildFields() self.sp.setFieldWithIndex(self.fields, i, 1.0) @@ -89,6 +93,9 @@ def buildDrivedEvolutionOperator(self): q = self.sp.fieldsAsStateVector(self.fields) A[:,i] = q[:] - self.fields = self.sp.buildFields() + self.fields = oldFields + + if reduceToEsentialDoF: + A = self.sp.reduceToEssentialDoF(A) return A \ No newline at end of file diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index a4e3232..94c4bcd 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -22,10 +22,10 @@ def __init__(self, mesh: Mesh1D): self.c0 = 1.0 self.tfsf = False self.source = None - + def TFSF_conditions(self, setup): - self.tfsf = True + self.tfsf = True self.source = setup["source"] self.left_TF_limit = (np.absolute(self.x - setup["left"])).argmin() self.right_TF_limit = (np.absolute(self.x - setup["right"])).argmin() @@ -35,7 +35,7 @@ def TFSF_conditions(self, setup): def buildFields(self): E = np.zeros(self.x.shape) H = np.zeros(self.xH.shape) - + if (self.source != None and self.tfsf): self.buildIncidentFields() @@ -46,16 +46,13 @@ def buildIncidentFields(self): self.Einc[:] = self.source(self.x[:]) self.Eprev = np.zeros(self.x.shape) - + self.Hinc = np.ndarray(self.xH.shape) self.Hinc[:] = self.source(self.xH[:] - 0.5*self.dt) - def get_minimum_node_distance(self): return np.min(self.dx) - - def computeRHSE(self, fields): H = fields['H'] E = fields['E'] @@ -66,18 +63,20 @@ def computeRHSE(self, fields): if self.tfsf == True: self.updateIncidentFieldE() - rhsE[self.left_TF_limit] += (1.0/self.dxH[0]) * self.Hinc[self.left_TF_limit-1] - rhsE[self.right_TF_limit] -= (1.0/self.dxH[0]) * self.Hinc[self.right_TF_limit ] + rhsE[self.left_TF_limit] += (1.0/self.dxH[0]) * \ + self.Hinc[self.left_TF_limit-1] + rhsE[self.right_TF_limit] -= (1.0/self.dxH[0]) * \ + self.Hinc[self.right_TF_limit] for bdr, label in self.mesh.boundary_label.items(): - + if bdr == "LEFT": if label == "PEC": rhsE[0] = 0.0 if label == "PMC": rhsE[0] = - (1.0/self.dxH[0]) * (2 * H[0]) - + if label == "Periodic": rhsE[0] = - (1.0/self.dxH[0]) * (H[0] - H[-1]) rhsE[-1] = rhsE[0] @@ -85,60 +84,46 @@ def computeRHSE(self, fields): if label == "Mur": rhsE[0] = E[1] + \ - (self.c0 * self.dt - self.dx[0]) / \ - (self.c0 * self.dt + self.dx[0]) * \ - (rhsE[1] - E[0]) - - rhsE[0] -= E[0] - rhsE[0] /= self.dt + (self.c0 * self.dt - self.dx[0]) / \ + (self.c0 * self.dt + self.dx[0]) * \ + (rhsE[1] - E[0]) - + rhsE[0] -= E[0] + rhsE[0] /= self.dt if bdr == "RIGHT": if label == "PEC": rhsE[-1] = 0.0 if label == "PMC": - rhsE[-1] = - (1.0/self.dxH[0]) * (-2 * H[-1]) - + rhsE[-1] = - (1.0/self.dxH[0]) * (-2 * H[-1]) + if label == "Periodic": rhsE[0] = - (1.0/self.dxH[0]) * (H[0] - H[-1]) rhsE[-1] = rhsE[0] - + if label == "Mur": rhsE[-1] = E[-2] + \ - (self.c0 * self.dt - self.dx[0]) / \ - (self.c0 * self.dt + self.dx[0]) * \ - (rhsE[-2] - E[-1]) - - rhsE[-1] -= E[-1] - rhsE[-1] /= self.dt - - - - # elif self.mesh.boundary_label == "PML": # [WIP] - # boundary_low = [0, 0] - # boundary_high = [0, 0] - - # rhsE[0] = boundary_low.pop(0) - # boundary_low.append(rhsE[1]) + (self.c0 * self.dt - self.dx[0]) / \ + (self.c0 * self.dt + self.dx[0]) * \ + (rhsE[-2] - E[-1]) - # rhsE[-1] = boundary_high.pop(0) - # boundary_high.append(rhsE[-2]) + rhsE[-1] -= E[-1] + rhsE[-1] /= self.dt return rhsE - - def computeRHSH(self, fields): E = fields['E'] rhsH = - (1.0/self.dx) * (E[1:] - E[:-1]) if self.tfsf == True: self.updateIncidentFieldH() - rhsH[self.left_TF_limit - 1] += (1.0/self.dx[0]) * self.Einc[self.left_TF_limit] - rhsH[self.right_TF_limit] -= (1.0/self.dx[0]) * self.Einc[self.right_TF_limit] + rhsH[self.left_TF_limit - + 1] += (1.0/self.dx[0]) * self.Einc[self.left_TF_limit] + rhsH[self.right_TF_limit] -= (1.0/self.dx[0]) * \ + self.Einc[self.right_TF_limit] return rhsH @@ -149,8 +134,9 @@ def computeRHS(self, fields): return {'E': rhsE, 'H': rhsH} def updateIncidentFieldE(self): - self.Einc[1:-1] = self.Einc[1:-1] - self.dt*(1.0/self.dxH) * (self.Hinc[1:] - self.Hinc[:-1]) - + self.Einc[1:-1] = self.Einc[1:-1] - self.dt * \ + (1.0/self.dxH) * (self.Hinc[1:] - self.Hinc[:-1]) + self.Einc[0] = \ self.Eprev[1] - \ (self.c0 * self.dt - self.dx[0]) / \ @@ -166,16 +152,15 @@ def updateIncidentFieldE(self): self.Eprev[:] = self.Einc[:] def updateIncidentFieldH(self): - self.Hinc = self.Hinc - self.dt*(1.0/self.dx) * (self.Einc[1:] - self.Einc[:-1]) - -#·································································································· + self.Hinc = self.Hinc - self.dt * \ + (1.0/self.dx) * (self.Einc[1:] - self.Einc[:-1]) def isStaggered(self): return True def number_of_nodes_per_element(self): - return 1 - + return 1 + def setFieldWithIndex(self, fields, i, val): NE = fields['E'].size if i < NE: @@ -184,70 +169,65 @@ def setFieldWithIndex(self, fields, i, val): fields['H'][i - NE] = val return fields - def buildEvolutionOperator(self): + def reduceToEssentialDoF(self, A): NE = self.buildFields()['E'].size - N = self.number_of_unknowns() - A = np.zeros((N, N)) - for i in range(N): - fields = self.buildFields() - self.setFieldWithIndex(fields, i, 1.0) - fieldsRHS = self.computeRHS(fields) - q0 = np.concatenate([fieldsRHS['E'], fieldsRHS['H']]) - A[:, i] = q0[:] - if self.mesh.boundary_label['LEFT'] == 'Periodic'\ - and self.mesh.boundary_label['RIGHT'] == 'Periodic' : + and self.mesh.boundary_label['RIGHT'] == 'Periodic': A = np.delete(A, NE-1, 0) - A[:,0] += A[:,NE-1] + A[:, 0] += A[:, NE-1] A = np.delete(A, NE-1, 1) elif self.mesh.boundary_label['LEFT'] == 'PEC'\ - and self.mesh.boundary_label['RIGHT'] == 'PEC': + and self.mesh.boundary_label['RIGHT'] == 'PEC': A = np.delete(A, NE-1, 0) A = np.delete(A, NE-1, 1) A = np.delete(A, 0, 0) - A = np.delete(A, 0, 1) + A = np.delete(A, 0, 1) else: - raise ValueError( "Periodic conditions must be ensured at both ends") + raise ValueError( + "Periodic conditions must be ensured at both ends") + + return A + + def buildEvolutionOperator(self, reduceToEssentialDoF=True): + N = self.number_of_unknowns() + A = np.zeros((N, N)) + for i in range(N): + fields = self.buildFields() + self.setFieldWithIndex(fields, i, 1.0) + fieldsRHS = self.computeRHS(fields) + q0 = np.concatenate([fieldsRHS['E'], fieldsRHS['H']]) + A[:, i] = q0[:] + + if reduceToEssentialDoF: + A = self.reduceToEssentialDoF(A) return A - def reorder_array(self, A, ordering): + def reorder_by_elements(self, A): # Assumes that the original array contains all DoF ordered as: # [ E_0, ..., E_{NE-1}, H_0, ..., H_{NH-1} ] N = A.shape[0] NE = len(self.x) NH = len(self.xH) if self.mesh.boundary_label['LEFT'] == 'Periodic' and \ - self.mesh.boundary_label['RIGHT'] == 'Periodic': + self.mesh.boundary_label['RIGHT'] == 'Periodic': NE -= 1 else: - raise ValueError( "Periodic conditions must be ensured at both ends") + raise ValueError( + "Periodic conditions must be ensured at both ends") if NE != NH: raise ValueError( "Unable to order by elements with different size fields.") N = NE + NH new_order = np.zeros(N, dtype=int) - 1 - if ordering == 'byElements': - for i in range(N): - if i < NE: - new_order[2*i] = i - else: - new_order[2*int(i - NE)+1] = i + + for i in range(N): + if i < NE: + new_order[2*i] = i + else: + new_order[2*int(i - NE)+1] = i + if (len(A.shape) == 1): A1 = [A[i] for i in new_order] elif (len(A.shape) == 2): A1 = [[A[i][j] for j in new_order] for i in new_order] return np.array(A1) - - def getEnergy(self, field): - raise ValueError("Not implemented") - Np = self.number_of_nodes_per_element() - K = self.mesh.number_of_elements() - assert field.shape == (Np, K) - energy = 0.0 - for k in range(K): - energy += np.inner( - field[:, k].dot(self.mass), - field[:, k]*self.jacobian[:, k] - ) - - return energy diff --git a/test/fd/test_fd1d.py b/test/fd/test_fd1d.py index 7782a18..7c8f605 100644 --- a/test/fd/test_fd1d.py +++ b/test/fd/test_fd1d.py @@ -63,7 +63,7 @@ def test_buildEvolutionOperator_sorting(): A = sp.buildEvolutionOperator() eigA, _ = np.linalg.eig(A) - A_by_elem = sp.reorder_array(A, 'byElements') + A_by_elem = sp.reorder_by_elements(A) eigA_by_elem, _ = np.linalg.eig(A_by_elem) assert A.shape == A_by_elem.shape diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 4a9555f..18f9493 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -939,7 +939,7 @@ def test_buildDrivedEvolutionOperator(): ) driver = MaxwellDriver(sp) - A = driver.buildDrivedEvolutionOperator() + A = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index dcb5810..c6494b3 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -27,7 +27,7 @@ def test_buildDrivedEvolutionOperator(): sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="PEC")) driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) - A = driver.buildDrivedEvolutionOperator() + A = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) @@ -41,6 +41,23 @@ def test_buildDrivedEvolutionOperator(): q = A.dot(q0) assert np.allclose(qExpected, q) + +def test_buildDrivedEvolutionOperator_reduced(): + K = 5 + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 5, boundary_label="Periodic")) + driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) + + A = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=True) + + # A = sp.reorder_by_elements(A) + # plt.matshow(A, cmap='RdGy') + # plt.colorbar(fraction=0.046, pad=0.04) + # for k in range(K): + # plt.vlines(k*2-0.5, -0.5, K*2-0.5, color='gray', linestyle='dashed') + # plt.hlines(k*2-0.5, -0.5, K*2-0.5, color='gray', linestyle='dashed') + # plt.show() + + assert np.allclose(np.abs(np.linalg.eig(A)[0]), 1.0) def test_fdtd_pec(): From ed857a1021f09b0696f9170ca6c7c135de18aed7 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 30 Apr 2024 10:42:14 +0200 Subject: [PATCH 02/25] changes to reoder_by_elements --- maxwell/dg/dg1d.py | 25 ++++++++----------------- test/dg/test_maxwell1d.py | 4 ++-- test/tests_fdtd/test_fdtd_1d.py | 17 ++++++++++------- 3 files changed, 20 insertions(+), 26 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index 44cb17c..eefaad7 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -215,29 +215,20 @@ def buildEvolutionOperator(self): A[:, i] = q0[:, 0] return A - def reorder_array(self, A, ordering): + def reorder_by_elements(self, A): # Assumes that the original array contains all DoF ordered as: # [ E_0, ..., E_{K-1}, H_0, ..., H_{K-1} ] N = A.shape[0] K = self.mesh.number_of_elements() Np = self.number_of_nodes_per_element() new_order = np.arange(N, dtype=int) - if ordering == 'byElements': - for i in range(N): - node = i % Np - elem = int(np.floor(i / Np)) % K - if i < N/2: - new_order[2*elem*Np+node] = i - else: - new_order[2*elem*Np+Np+node] = i - if ordering == 'interleaved': - for i in range(N): - node = i % Np - elem = int(np.floor(i / Np)) % K - if i < N/2: - new_order[2*elem*Np+node*2] = i - else: - new_order[2*elem*Np+node*2+1] = i + for i in range(N): + node = i % Np + elem = int(np.floor(i / Np)) % K + if i < N/2: + new_order[2*elem*Np+node] = i + else: + new_order[2*elem*Np+Np+node] = i if (len(A.shape) == 1): A1 = [A[i] for i in new_order] elif (len(A.shape) == 2): diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index 13d8382..879c8f4 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -22,7 +22,7 @@ def test_buildEvolutionOperator_PEC(): m = Mesh1D(0, 1, 5, boundary_label='PEC') sp = DG1D(1, m, "Centered") A = sp.buildEvolutionOperator() - A = sp.reorder_array(A, 'byElements') + A = sp.reorder_by_elements(A) M = sp.buildGlobalMassMatrix() @@ -51,7 +51,7 @@ def test_buildEvolutionOperator_sorting(): A = sp.buildEvolutionOperator() eigA, _ = np.linalg.eig(A) - A_by_elem = sp.reorder_array(A, 'byElements') + A_by_elem = sp.reorder_by_elements(A) eigA_by_elem, _ = np.linalg.eig(A_by_elem) assert A.shape == A_by_elem.shape diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index c6494b3..4a36416 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -17,7 +17,7 @@ def plot(sp, driver): plt.ylim(-1, 1) plt.title(driver.timeIntegrator.time) plt.grid(which='both') - plt.pause(0.01) + plt.pause(0.1) plt.cla() # ······················································ @@ -44,10 +44,11 @@ def test_buildDrivedEvolutionOperator(): def test_buildDrivedEvolutionOperator_reduced(): K = 5 - sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 5, boundary_label="Periodic")) - driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) - - A = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=True) + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="PEC")) + + A_0_9 = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.9).buildDrivedEvolutionOperator() + A_1_0 = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0).buildDrivedEvolutionOperator() + A_1_01 = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.01).buildDrivedEvolutionOperator() # A = sp.reorder_by_elements(A) # plt.matshow(A, cmap='RdGy') @@ -57,7 +58,9 @@ def test_buildDrivedEvolutionOperator_reduced(): # plt.hlines(k*2-0.5, -0.5, K*2-0.5, color='gray', linestyle='dashed') # plt.show() - assert np.allclose(np.abs(np.linalg.eig(A)[0]), 1.0) + assert np.allclose(np.abs(np.linalg.eig(A_0_9)[0]), 1.0) + assert np.allclose(np.abs(np.linalg.eig(A_1_0)[0]), 1.0) + assert np.any(np.abs(np.linalg.eig(A_1_01)[0]) - 1.0 > 0) def test_fdtd_pec(): @@ -78,7 +81,7 @@ def test_fdtd_pec(): def test_fdtd_periodic(): - sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="Periodic")) + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 5, boundary_label="Periodic")) driver = MaxwellDriver(sp, timeIntegratorType='LF2') s0 = 0.25 From 29c4c783fd5c625fbac375fbcabbc1e2b96ec2e9 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 1 May 2024 13:17:44 +0200 Subject: [PATCH 03/25] Fixes Verlet tests --- maxwell/integrators/LF2V.py | 16 +------ test/dg/test_maxwell1d.py | 17 ++++++++ test/test_dgtd_1d.py | 85 ++++++++++++++++++++++++------------- 3 files changed, 75 insertions(+), 43 deletions(-) diff --git a/maxwell/integrators/LF2V.py b/maxwell/integrators/LF2V.py index 834bc76..867beaf 100644 --- a/maxwell/integrators/LF2V.py +++ b/maxwell/integrators/LF2V.py @@ -12,22 +12,10 @@ def step(self, fields, dt): E = fields['E'] H = fields['H'] - # #Velocity Verlet - # self.time += dt - # H += 0.5*dt*self.sp.computeRHSH(fields) - - # self.time += dt - # E += dt*self.sp.computeRHSE(fields) - # H += 0.5*dt*self.sp.computeRHSH(fields) - - - #Position Verlet - E += 0.5*dt*self.sp.computeRHSE(fields) - - self.time += dt - H += 0.5*dt*self.sp.computeRHSH(fields) + H += dt*self.sp.computeRHSH(fields) E += 0.5*dt*self.sp.computeRHSE(fields) + self.time += dt # #Verlet Algorithm diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index 879c8f4..564802d 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -18,6 +18,23 @@ def test_get_energy_N1(): assert np.isclose(sp.getEnergy(fields['E']), 1.0, rtol=1e-9) +def test_energy_with_operators(): + m = Mesh1D(0, 1, 10) + sp = DG1D(1, m) + + fields = sp.buildFields() + fields['E'].fill(1.0) + fields['H'].fill(2.0) + + expectedEnergy = sp.getEnergy(fields['E']) + sp.getEnergy(fields['H']) + + q = sp.fieldsAsStateVector(fields) + Mg = sp.buildGlobalMassMatrix() + energy = q.T.dot(Mg).dot(q) + + assert np.isclose(expectedEnergy, energy) + + def test_buildEvolutionOperator_PEC(): m = Mesh1D(0, 1, 5, boundary_label='PEC') sp = DG1D(1, m, "Centered") diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 18f9493..95b996d 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -215,17 +215,16 @@ def test_pec_centered_lf2(): # plt.cla() -@pytest.mark.skip(reason="Doesn't work. Deactivated to pass automated tests.") def test_pec_centered_lf2v(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), fluxType="Centered" ) - driver = MaxwellDriver(sp, timeIntegratorType='LF2V') + driver = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.8) final_time = 1.999 - s0 = 0.15 + s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] @@ -235,18 +234,46 @@ def test_pec_centered_lf2v(): R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), -finalFieldE.reshape(1, finalFieldE.size)) - # assert R[0,1] > 0.9999 + assert R[0,1] > 0.9999 - driver['E'][:] = initialFieldE[:] - for _ in range(1000): - driver.step() - plt.plot(sp.x, driver['E'], 'b') - plt.ylim(-1, 1) - plt.grid(which='both') - plt.pause(0.01) - plt.cla() + # driver['E'][:] = initialFieldE[:] + # for _ in range(1000): + # driver.step() + # plt.plot(sp.x, driver['E'], 'b') + # plt.plot(sp.x, driver['H'], 'r') + # plt.ylim(-1, 1) + # plt.grid(which='both') + # plt.pause(0.01) + # plt.cla() +def test_pec_centered_lf2_and_lf2v_are_equivalent(): + sp = DG1D( + n_order=3, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), + fluxType="Centered" + ) + drLF = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.8) + drVe = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.8) + + s0 = 0.25 + initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) + drLF['E'][:] = initialFieldE + drVe['E'][:] = initialFieldE + + for _ in range(200): + drLF.step() + drVe.step() + assert np.allclose(drLF['H'] - drVe['H'], 0.0) + + # plt.plot(sp.x, drLF['H'], 'b') + # plt.plot(sp.x, drVe['H'], 'g') + # plt.ylim(-1, 1) + # plt.grid(which='both') + # plt.pause(0.01) + # plt.cla() + + def test_pec_centered_ibe(): sp = DG1D( n_order=3, @@ -395,10 +422,10 @@ def test_energy_evolution_centered(): totalEnergy = energyE + energyH assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 3e-5 - # plt.plot(energyE, 'b') - # plt.plot(energyH, 'r') - # plt.plot(totalEnergy, 'g') - # plt.show() + plt.plot(energyE, 'b') + plt.plot(energyH, 'r') + plt.plot(totalEnergy, 'g') + plt.show() def test_energy_evolution_centered_lf2(): @@ -435,7 +462,6 @@ def test_energy_evolution_centered_lf2(): # plt.show() -@pytest.mark.skip(reason="Doesn't work. Deactivated to pass automated tests.") def test_energy_evolution_centered_lf2v(): sp = DG1D( n_order=2, @@ -443,7 +469,7 @@ def test_energy_evolution_centered_lf2v(): fluxType="Centered" ) - driver = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.85) + driver = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.7) driver['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) Nsteps = 500 @@ -452,22 +478,23 @@ def test_energy_evolution_centered_lf2v(): for n in range(Nsteps): energyE[n] = sp.getEnergy(driver['E']) energyH[n] = sp.getEnergy(driver['H']) - # plt.plot(sp.x, driver['E'], 'b') - # plt.ylim(-1,1) - # plt.grid(which='both') - # plt.pause(0.1) - # plt.cla() + plt.plot(sp.x, driver['E'],'b') + plt.plot(sp.x, driver['H'],'r') + plt.ylim(-1,1) + plt.grid(which='both') + plt.pause(0.1) + plt.cla() driver.step() totalEnergy = energyE + energyH assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 3e-3 - # plt.figure() - # # plt.plot(energyE) - # # plt.plot(energyH) - # plt.plot(totalEnergy) - # plt.grid() - # plt.show() + plt.figure() + plt.plot(energyE) + plt.plot(energyH) + plt.plot(totalEnergy) + plt.grid() + plt.show() def test_periodic_tested(): From 627047d5a3f18831d318f25ec07a2b0d2b6a9eb2 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 1 May 2024 14:04:09 +0200 Subject: [PATCH 04/25] Experimenting with operators --- maxwell/driver.py | 2 +- maxwell/integrators/LF2V.py | 15 ------ test/test_dgtd_1d.py | 94 ++++++++++++++++++++++++++++--------- 3 files changed, 72 insertions(+), 39 deletions(-) diff --git a/maxwell/driver.py b/maxwell/driver.py index 26deae4..fc44a93 100644 --- a/maxwell/driver.py +++ b/maxwell/driver.py @@ -95,7 +95,7 @@ def buildDrivedEvolutionOperator(self, reduceToEsentialDoF=True): self.fields = oldFields - if reduceToEsentialDoF: + if reduceToEsentialDoF and self.sp.isStaggered(): A = self.sp.reduceToEssentialDoF(A) return A \ No newline at end of file diff --git a/maxwell/integrators/LF2V.py b/maxwell/integrators/LF2V.py index 867beaf..4acfff0 100644 --- a/maxwell/integrators/LF2V.py +++ b/maxwell/integrators/LF2V.py @@ -18,21 +18,6 @@ def step(self, fields, dt): self.time += dt - # #Verlet Algorithm - # self.time += dt - # E += dt*self.sp.computeRHSE(fields) - # H += 0.5*dt*self.sp.computeRHSH(fields) - # E_old = E - - # self.time += dt/2 - - # H += 0.5*dt*self.sp.computeRHSH(fields) - # E += dt*self.sp.computeRHSE(fields) - # E_new = E - - # E += 2*E_new - E_old + dt**2*self.sp.computeRHSH(fields) - - diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 95b996d..3474b5b 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -403,7 +403,7 @@ def test_energy_evolution_centered(): dissipate because of the LSERK4 time integration. ''' sp = DG1D( - n_order=5, + n_order=2, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), fluxType="Centered" ) @@ -420,13 +420,38 @@ def test_energy_evolution_centered(): driver.step() totalEnergy = energyE + energyH - assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 3e-5 + # Expect some dissipation due to LSERK4. + assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 - plt.plot(energyE, 'b') - plt.plot(energyH, 'r') - plt.plot(totalEnergy, 'g') - plt.show() + # plt.plot(energyE, 'b') + # plt.plot(energyH, 'r') + # plt.plot(totalEnergy, 'g') + # plt.show() +def test_energy_evolution_with_operators(): + sp = DG1D( + n_order=2, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), + fluxType="Centered" + ) + dr = MaxwellDriver(sp, timeIntegratorType='LSERK4', CFL=0.7) + dr['H'][:] = np.exp(-sp.x**2/(2*0.25**2)) + + q0 = sp.fieldsAsStateVector(dr.fields) + G = dr.buildDrivedEvolutionOperator() + Mg = sp.buildGlobalMassMatrix() + q = G.dot(q0) + + initialEnergy = q0.T.dot(Mg).dot(q0) + finalEnergy = q.T.dot(Mg).dot(q) + + expectedPower = (finalEnergy - initialEnergy)/dr.dt + + power = (1/dr.dt) * q0.T.dot(G.T.dot(Mg).dot(G) - Mg).dot(q0) + + assert np.isclose(power, expectedPower) + + def test_energy_evolution_centered_lf2(): sp = DG1D( @@ -454,7 +479,6 @@ def test_energy_evolution_centered_lf2(): totalEnergy = energyE + energyH assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 0.01 - # plt.figure() # plt.plot(energyE) # plt.plot(energyH) # plt.plot(totalEnergy) @@ -464,7 +488,7 @@ def test_energy_evolution_centered_lf2(): def test_energy_evolution_centered_lf2v(): sp = DG1D( - n_order=2, + n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), fluxType="Centered" ) @@ -478,23 +502,23 @@ def test_energy_evolution_centered_lf2v(): for n in range(Nsteps): energyE[n] = sp.getEnergy(driver['E']) energyH[n] = sp.getEnergy(driver['H']) - plt.plot(sp.x, driver['E'],'b') - plt.plot(sp.x, driver['H'],'r') - plt.ylim(-1,1) - plt.grid(which='both') - plt.pause(0.1) - plt.cla() + # plt.plot(sp.x, driver['E'],'b') + # plt.plot(sp.x, driver['H'],'r') + # plt.ylim(-1,1) + # plt.grid(which='both') + # plt.pause(0.1) + # plt.cla() driver.step() totalEnergy = energyE + energyH - assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 3e-3 + # assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 3e-3 - plt.figure() - plt.plot(energyE) - plt.plot(energyH) - plt.plot(totalEnergy) - plt.grid() - plt.show() + # plt.figure() + # plt.plot(energyE) + # plt.plot(energyH) + # plt.plot(totalEnergy) + # plt.grid() + # plt.show() def test_periodic_tested(): @@ -966,7 +990,31 @@ def test_buildDrivedEvolutionOperator(): ) driver = MaxwellDriver(sp) - A = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) + G = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) + + s0 = 0.25 + initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) + driver['E'][:] = initialFieldE[:] + + q0 = sp.fieldsAsStateVector(driver.fields) + + driver.step() + qExpected = sp.fieldsAsStateVector(driver.fields) + + q = G.dot(q0) + + assert np.allclose(qExpected, q) + + +def test_buildDrivedEvolutionOperator_LF2V(): + sp = DG1D( + n_order=2, + mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), + fluxType="Centered" + ) + driver = MaxwellDriver(sp, timeIntegratorType='LF2V') + + G = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) @@ -977,6 +1025,6 @@ def test_buildDrivedEvolutionOperator(): driver.step() qExpected = sp.fieldsAsStateVector(driver.fields) - q = A.dot(q0) + q = G.dot(q0) assert np.allclose(qExpected, q) \ No newline at end of file From dcf8419a9edd32026a9487305c4a6823dafd0895 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 1 May 2024 19:46:16 +0200 Subject: [PATCH 05/25] Adds energy operators test --- test/test_dgtd_1d.py | 91 ++++++++++++++++---------------------------- 1 file changed, 32 insertions(+), 59 deletions(-) diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 3474b5b..2c814cc 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -23,6 +23,7 @@ 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(): sp = DG1D( n_order=5, @@ -234,7 +235,7 @@ def test_pec_centered_lf2v(): R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), -finalFieldE.reshape(1, finalFieldE.size)) - assert R[0,1] > 0.9999 + assert R[0, 1] > 0.9999 # driver['E'][:] = initialFieldE[:] # for _ in range(1000): @@ -255,24 +256,24 @@ def test_pec_centered_lf2_and_lf2v_are_equivalent(): ) drLF = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.8) drVe = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.8) - + s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) drLF['E'][:] = initialFieldE drVe['E'][:] = initialFieldE - + for _ in range(200): drLF.step() drVe.step() - assert np.allclose(drLF['H'] - drVe['H'], 0.0) - + assert np.allclose(drLF['H'] - drVe['H'], 0.0) + # plt.plot(sp.x, drLF['H'], 'b') # plt.plot(sp.x, drVe['H'], 'g') # plt.ylim(-1, 1) # plt.grid(which='both') # plt.pause(0.01) # plt.cla() - + def test_pec_centered_ibe(): sp = DG1D( @@ -367,35 +368,6 @@ def test_periodic_centered_dirk2(): # plt.pause(0.001) # plt.cla() -# def test_pec_centered_iglrk4(): -# sp = Maxwell1D( -# n_order = 3, -# mesh = Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), -# fluxType="Upwind" -# ) -# driver = MaxwellDriver(sp, timeIntegratorType='IGLRK4', CFL=2) - -# final_time = 3.999 -# s0 = 0.25 -# initialField = np.exp(-(sp.x)**2/(2*s0**2)) - -# driver['E'][:] = initialField[:] -# driver['H'][:] = initialField[:] -# finalFieldE = driver['E'] - -# driver.run_until(final_time) - -# driver['E'][:] = initialField[:] -# driver['H'][:] = initialField[:] -# for _ in range(500): -# driver.step() -# plt.plot(sp.x, driver['E'],'b') -# plt.plot(sp.x, driver['H'],'r') -# plt.ylim(-1, 1) -# plt.grid(which='both') -# plt.pause(0.001) -# plt.cla() - def test_energy_evolution_centered(): ''' @@ -421,37 +393,38 @@ def test_energy_evolution_centered(): totalEnergy = energyE + energyH # Expect some dissipation due to LSERK4. - assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 + assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 # plt.plot(energyE, 'b') # plt.plot(energyH, 'r') # plt.plot(totalEnergy, 'g') # plt.show() + def test_energy_evolution_with_operators(): sp = DG1D( n_order=2, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), + mesh=Mesh1D(0, 1.0, 10, boundary_label="Periodic"), fluxType="Centered" ) - dr = MaxwellDriver(sp, timeIntegratorType='LSERK4', CFL=0.7) - dr['H'][:] = np.exp(-sp.x**2/(2*0.25**2)) - + dr = MaxwellDriver(sp) + dr['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) + q0 = sp.fieldsAsStateVector(dr.fields) G = dr.buildDrivedEvolutionOperator() Mg = sp.buildGlobalMassMatrix() q = G.dot(q0) - + initialEnergy = q0.T.dot(Mg).dot(q0) finalEnergy = q.T.dot(Mg).dot(q) - + expectedPower = (finalEnergy - initialEnergy)/dr.dt + PG = (1/dr.dt)*(G.T.dot(Mg).dot(G) - Mg) - power = (1/dr.dt) * q0.T.dot(G.T.dot(Mg).dot(G) - Mg).dot(q0) - + power = q0.T.dot(PG).dot(q0) + assert np.isclose(power, expectedPower) - - + def test_energy_evolution_centered_lf2(): sp = DG1D( @@ -989,23 +962,23 @@ def test_buildDrivedEvolutionOperator(): fluxType="Centered" ) driver = MaxwellDriver(sp) - + G = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) - + s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] - + q0 = sp.fieldsAsStateVector(driver.fields) - + driver.step() qExpected = sp.fieldsAsStateVector(driver.fields) - + q = G.dot(q0) assert np.allclose(qExpected, q) - - + + def test_buildDrivedEvolutionOperator_LF2V(): sp = DG1D( n_order=2, @@ -1013,18 +986,18 @@ def test_buildDrivedEvolutionOperator_LF2V(): fluxType="Centered" ) driver = MaxwellDriver(sp, timeIntegratorType='LF2V') - + G = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) - + s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] - + q0 = sp.fieldsAsStateVector(driver.fields) - + driver.step() qExpected = sp.fieldsAsStateVector(driver.fields) - + q = G.dot(q0) - assert np.allclose(qExpected, q) \ No newline at end of file + assert np.allclose(qExpected, q) From 5267d247d2184fd0cbefb13b32af0406c549d42e Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 1 May 2024 19:49:27 +0200 Subject: [PATCH 06/25] Adds buildPowerOperator --- maxwell/driver.py | 7 ++++++- test/test_dgtd_1d.py | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/maxwell/driver.py b/maxwell/driver.py index fc44a93..12330a7 100644 --- a/maxwell/driver.py +++ b/maxwell/driver.py @@ -98,4 +98,9 @@ def buildDrivedEvolutionOperator(self, reduceToEsentialDoF=True): if reduceToEsentialDoF and self.sp.isStaggered(): A = self.sp.reduceToEssentialDoF(A) - return A \ No newline at end of file + return A + + def buildPowerOperator(self): + G = self.buildDrivedEvolutionOperator() + Mg = self.sp.buildGlobalMassMatrix() + return (1/self.dt)*(G.T.dot(Mg).dot(G) - Mg) \ No newline at end of file diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 2c814cc..d1413a6 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -419,8 +419,8 @@ def test_energy_evolution_with_operators(): finalEnergy = q.T.dot(Mg).dot(q) expectedPower = (finalEnergy - initialEnergy)/dr.dt - PG = (1/dr.dt)*(G.T.dot(Mg).dot(G) - Mg) + PG = dr.buildPowerOperator() power = q0.T.dot(PG).dot(q0) assert np.isclose(power, expectedPower) From 7b431d5ecd7b378e583b0c25ae6f1ceffab1842e Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 7 May 2024 09:33:33 +0200 Subject: [PATCH 07/25] Reorganizing operatos --- maxwell/driver.py | 22 +++++++- maxwell/integrators/LSERK4.py | 8 ++- maxwell/spatialDiscretization.py | 29 ++++++++--- test/test_dgtd_1d.py | 89 ++++++++++++++++++++------------ 4 files changed, 102 insertions(+), 46 deletions(-) diff --git a/maxwell/driver.py b/maxwell/driver.py index 12330a7..8da147b 100644 --- a/maxwell/driver.py +++ b/maxwell/driver.py @@ -103,4 +103,24 @@ def buildDrivedEvolutionOperator(self, reduceToEsentialDoF=True): def buildPowerOperator(self): G = self.buildDrivedEvolutionOperator() Mg = self.sp.buildGlobalMassMatrix() - return (1/self.dt)*(G.T.dot(Mg).dot(G) - Mg) \ No newline at end of file + return (1/self.dt)*(G.T.dot(Mg).dot(G) - Mg) + + + def buildCausallyConnectedOperators(self, neighbors=-1): + if neighbors == -1: + neighs = self.timeIntegrator.N_STAGES + else: + neighs = neighbors + + local_indices, neigh_indices = self.sp.buildLocalAndNeighborIndices(neighs) + + G = self.sp.reorder_by_elements(self.buildDrivedEvolutionOperator()) + Mg = self.sp.buildGlobalMassMatrix() + A = G[local_indices][:,local_indices] + B = G[local_indices][:,neigh_indices] + C = G[neigh_indices][:,local_indices] + D = G[neigh_indices][:,neigh_indices] + Mk = Mg[local_indices][:,local_indices] + Mn = Mg[neigh_indices][:,neigh_indices] + + return A, B, C, D, Mk, Mn \ No newline at end of file diff --git a/maxwell/integrators/LSERK4.py b/maxwell/integrators/LSERK4.py index 76d13a6..3fd59f5 100644 --- a/maxwell/integrators/LSERK4.py +++ b/maxwell/integrators/LSERK4.py @@ -26,10 +26,8 @@ def __init__(self, sp: SpatialDiscretization, fields): def step(self, fields, dt): for s in range(0, self.N_STAGES): fieldsRHS = self.sp.computeRHS(fields) - for l, f in fieldsRHS.items(): - self.fieldsRes[l] = self.A[s]*self.fieldsRes[l] + dt*f - - for l, f in fields.items(): - fields[l] += self.B[s] * self.fieldsRes[l] + for label, f in fieldsRHS.items(): + self.fieldsRes[label] = self.A[s]*self.fieldsRes[label] + dt*f + fields[label] += self.B[s] * self.fieldsRes[label] self.time += dt diff --git a/maxwell/spatialDiscretization.py b/maxwell/spatialDiscretization.py index b4f850e..b7f3082 100644 --- a/maxwell/spatialDiscretization.py +++ b/maxwell/spatialDiscretization.py @@ -1,38 +1,51 @@ import numpy as np + class SpatialDiscretization(): def __init__(self, mesh): self.mesh = mesh return - + def get_mesh(self): return self.mesh - + def isStaggered(self): return False - + def dimension(self): return 1 def fieldsAsStateVector(self, fields): q = np.array([]) for f in fields.values(): - q = np.append(q, f.reshape(f.size,1, order='F')) + q = np.append(q, f.reshape(f.size, 1, order='F')) return q - + def buildStateVector(self): fields = self.buildFields() q0 = np.array([]) for f in fields.values(): q0 = np.append(q0, f) return q0 - + def buildImpulseStateVector(self, i): q = self.buildStateVector() q[i] = 1.0 return q - + + def buildLocalAndNeighborIndices(self, neighs): + Np = self.number_of_nodes_per_element() + N = self.number_of_unknowns() + + local_indices = np.arange(0, 2*Np) + right_neigh_indices = local_indices[-1] + 1 + np.arange(0, 2*Np*neighs) + left_neigh_indices = np.sort( + local_indices[0] - np.arange(0, 2*Np*neighs) - 1) % N + neigh_indices = np.sort(np.concatenate( + (left_neigh_indices, right_neigh_indices))) + + return local_indices, neigh_indices + def number_of_unknowns(self): return len(self.buildStateVector()) - diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index d1413a6..c4884cd 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -275,6 +275,7 @@ def test_pec_centered_lf2_and_lf2v_are_equivalent(): # plt.cla() + def test_pec_centered_ibe(): sp = DG1D( n_order=3, @@ -368,6 +369,25 @@ def test_periodic_centered_dirk2(): # plt.pause(0.001) # plt.cla() +def test_buildCausallyConnectedOperators(): + sp = DG1D( + n_order=1, + mesh=Mesh1D(-1.0, 1.0, 11, boundary_label="Periodic"), + fluxType="Upwind" + ) + dr = MaxwellDriver(sp) + + G = sp.reorder_by_elements(dr.buildDrivedEvolutionOperator()) + A,B,C,D,Mk,Mn = dr.buildCausallyConnectedOperators() + + G1 = np.concatenate([A,B], axis=1) + G2 = np.concatenate([C,D], axis=1) + G_reassembled = np.concatenate([G1, G2]) + + assert np.all(G == G_reassembled) + + + def test_energy_evolution_centered(): ''' @@ -375,9 +395,9 @@ def test_energy_evolution_centered(): dissipate because of the LSERK4 time integration. ''' sp = DG1D( - n_order=2, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + n_order=1, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), + fluxType="Upwind" ) driver = MaxwellDriver(sp) @@ -393,22 +413,24 @@ def test_energy_evolution_centered(): totalEnergy = energyE + energyH # Expect some dissipation due to LSERK4. - assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 + # assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 - # plt.plot(energyE, 'b') - # plt.plot(energyH, 'r') - # plt.plot(totalEnergy, 'g') + # plt.plot(energyE, '.-b') + # plt.plot(energyH, '.-r') + # plt.plot(totalEnergy, '.-g') # plt.show() def test_energy_evolution_with_operators(): sp = DG1D( - n_order=2, + n_order=1, mesh=Mesh1D(0, 1.0, 10, boundary_label="Periodic"), - fluxType="Centered" + fluxType="Upwind" + ) dr = MaxwellDriver(sp) dr['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) + dr['H'][:] = np.exp(-sp.x**2/(2*0.25**2)) q0 = sp.fieldsAsStateVector(dr.fields) G = dr.buildDrivedEvolutionOperator() @@ -428,33 +450,35 @@ def test_energy_evolution_with_operators(): def test_energy_evolution_centered_lf2(): sp = DG1D( - n_order=2, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), + n_order=3, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), fluxType="Centered" ) - driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.7) - driver['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) - - Nsteps = 1500 + dr = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.7) + dr['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) + dr['H'][:] = np.exp(-(sp.x-dr.dt/2)**2/(2*0.25**2)) + Nsteps = 500 energyE = np.zeros(Nsteps) energyH = np.zeros(Nsteps) for n in range(Nsteps): - energyE[n] = sp.getEnergy(driver['E']) - energyH[n] = sp.getEnergy(driver['H']) - # plt.plot(sp.x, driver['E'], 'b') + energyE[n] = sp.getEnergy(dr['E']) + energyH[n] = sp.getEnergy(dr['H']) + # plt.plot(sp.x, dr['E'], 'b') + # plt.plot(sp.x, dr['H'], 'r') # plt.ylim(-1,1) # plt.grid(which='both') # plt.pause(0.1) # plt.cla() - driver.step() + dr.step() - totalEnergy = energyE + energyH - assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 0.01 + totalEnergy = energyE[1:] + (energyH[:-1]+ energyH[1:])*.5 + + assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 1e-6 # plt.plot(energyE) # plt.plot(energyH) - # plt.plot(totalEnergy) + # plt.plot(totalEnergy, '-b') # plt.grid() # plt.show() @@ -462,29 +486,30 @@ def test_energy_evolution_centered_lf2(): def test_energy_evolution_centered_lf2v(): sp = DG1D( n_order=3, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), fluxType="Centered" ) - driver = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.7) - driver['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) + dr = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.7) + dr['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) + dr['H'][:] = np.exp(-sp.x**2/(2*0.25**2)) Nsteps = 500 energyE = np.zeros(Nsteps) energyH = np.zeros(Nsteps) for n in range(Nsteps): - energyE[n] = sp.getEnergy(driver['E']) - energyH[n] = sp.getEnergy(driver['H']) - # plt.plot(sp.x, driver['E'],'b') - # plt.plot(sp.x, driver['H'],'r') + energyE[n] = sp.getEnergy(dr['E']) + energyH[n] = sp.getEnergy(dr['H']) + # plt.plot(sp.x, dr['E'],'b') + # plt.plot(sp.x, dr['H'],'r') # plt.ylim(-1,1) # plt.grid(which='both') # plt.pause(0.1) # plt.cla() - driver.step() + dr.step() totalEnergy = energyE + energyH - # assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 3e-3 + assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 1e-6 # plt.figure() # plt.plot(energyE) @@ -959,7 +984,7 @@ def test_buildDrivedEvolutionOperator(): sp = DG1D( n_order=2, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxType="Upwind" ) driver = MaxwellDriver(sp) From 79637cb721bcd95e8e0287fadd435c7141a6c9f9 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 7 May 2024 11:23:29 +0200 Subject: [PATCH 08/25] Adds test for causally connected operator --- maxwell/driver.py | 4 ++-- maxwell/spatialDiscretization.py | 6 +++--- test/test_dgtd_1d.py | 30 ++++++++++++++++++++++++++++++ 3 files changed, 35 insertions(+), 5 deletions(-) diff --git a/maxwell/driver.py b/maxwell/driver.py index 8da147b..96985fd 100644 --- a/maxwell/driver.py +++ b/maxwell/driver.py @@ -106,13 +106,13 @@ def buildPowerOperator(self): return (1/self.dt)*(G.T.dot(Mg).dot(G) - Mg) - def buildCausallyConnectedOperators(self, neighbors=-1): + def buildCausallyConnectedOperators(self, element=0, neighbors=-1): if neighbors == -1: neighs = self.timeIntegrator.N_STAGES else: neighs = neighbors - local_indices, neigh_indices = self.sp.buildLocalAndNeighborIndices(neighs) + local_indices, neigh_indices = self.sp.buildLocalAndNeighborIndices(element, neighs) G = self.sp.reorder_by_elements(self.buildDrivedEvolutionOperator()) Mg = self.sp.buildGlobalMassMatrix() diff --git a/maxwell/spatialDiscretization.py b/maxwell/spatialDiscretization.py index b7f3082..1b16227 100644 --- a/maxwell/spatialDiscretization.py +++ b/maxwell/spatialDiscretization.py @@ -34,11 +34,11 @@ def buildImpulseStateVector(self, i): q[i] = 1.0 return q - def buildLocalAndNeighborIndices(self, neighs): + def buildLocalAndNeighborIndices(self, element, neighs): Np = self.number_of_nodes_per_element() N = self.number_of_unknowns() - - local_indices = np.arange(0, 2*Np) + k = element + local_indices = np.arange(k*2*Np, (k+1)*2*Np) right_neigh_indices = local_indices[-1] + 1 + np.arange(0, 2*Np*neighs) left_neigh_indices = np.sort( local_indices[0] - np.arange(0, 2*Np*neighs) - 1) % N diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index c4884cd..16df185 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -385,9 +385,39 @@ def test_buildCausallyConnectedOperators(): G_reassembled = np.concatenate([G1, G2]) assert np.all(G == G_reassembled) + + +def test_buildCausallyConnectedOperators_2(): + sp = DG1D( + n_order=1, + mesh=Mesh1D(-1.0, 1.0, 15, boundary_label="Periodic"), + fluxType="Upwind" + ) + dr = MaxwellDriver(sp) + G = sp.reorder_by_elements(dr.buildDrivedEvolutionOperator()) + A,B,C,D,Mk,Mn = dr.buildCausallyConnectedOperators(element=5) + G1 = np.concatenate([A, B], axis=1) + G2 = np.concatenate([C, D], axis=1) + G_reassembled = np.concatenate([G1, G2]) + Nk = A.shape[0] + Nn = int(B.shape[1]/2) + new_order = np.concatenate([np.arange(Nk, Nk+Nn), np.arange(0, Nk), np.arange(Nn+Nk, 2*Nn+Nk)]) + G_f = np.array([[G_reassembled[i][j] for j in new_order] for i in new_order]) + assert np.all(G[:G_f.shape[0], :G_f.shape[1]] == G_f) +def test_local_causal_operator(): + sp = DG1D( + n_order=1, + mesh=Mesh1D(-1.0, 1.0, 11, boundary_label="Periodic"), + fluxType="Centered" + ) + dr = MaxwellDriver(sp) + G = sp.reorder_by_elements(dr.buildDrivedEvolutionOperator()) + L = sp.reorder + + def test_energy_evolution_centered(): ''' From 8d82060af63633f82f346cd21ed4e141a4447ac2 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 7 May 2024 20:45:05 +0200 Subject: [PATCH 09/25] [WIP] changing upwind/centered flux to penalty. --- maxwell/dg/dg1d.py | 42 +++----------- maxwell/integrators/LSERK4.py | 1 + test/test_dgtd_1d.py | 105 ++++++++++++++++++++++------------ 3 files changed, 78 insertions(+), 70 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index eefaad7..e5f77fd 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -6,20 +6,16 @@ class DG1D(SpatialDiscretization): - def __init__(self, n_order: int, mesh: Mesh1D, fluxType="Upwind"): + def __init__(self, n_order: int, mesh: Mesh1D, fluxPenalty=1.0): SpatialDiscretization.__init__(self, mesh) assert n_order > 0 self.n_order = n_order - assert fluxType == "Upwind" or fluxType == "Centered" - self.fluxType = fluxType - + self.alpha = fluxPenalty self.n_faces = 2 self.n_fp = 1 - alpha = 0 - beta = 0 # Set up material parameters self.epsilon = np.ones(mesh.number_of_elements()) @@ -32,7 +28,7 @@ def __init__(self, n_order: int, mesh: Mesh1D, fluxType="Upwind"): self.vmap_m, self.vmap_p, self.vmap_b, self.map_b = build_maps( n_order, self.x, etoe, etof) - r = jacobiGL(alpha, beta, n_order) + r = jacobiGL(0, 0, n_order) self.fmask, self.fmask_1, self.fmask_2 = buildFMask(r) self.mass = mass_matrix(n_order, r) @@ -92,6 +88,9 @@ def fieldsOnBoundaryConditions(self, E, H): elif label == "PMC": Hbc = - H.transpose().take(self.vmap_b) Ebc = E.transpose().take(self.vmap_b) + elif label == "NULL": + Hbc = - H.transpose().take(self.vmap_b) + Ebc = - E.transpose().take(self.vmap_b) elif label == "SMA": Hbc = H.transpose().take(self.vmap_b) * 0.0 Ebc = E.transpose().take(self.vmap_b) * 0.0 @@ -104,39 +103,14 @@ def fieldsOnBoundaryConditions(self, E, H): def computeFluxE(self, E, H): dE, dH = self.computeJumps(E, H) - - if self.fluxType == "Upwind": - flux_E = 1/self.Z_imp_sum*(self.nx*self.Z_imp_p*dH-dE) - elif self.fluxType == "Centered": - flux_E = 1/self.Z_imp_sum*(self.nx*self.Z_imp_p*dH) - else: - raise ValueError("Invalid fluxType label") + flux_E = 1/self.Z_imp_sum*(self.nx*self.Z_imp_p*dH-self.alpha*dE) return flux_E def computeFluxH(self, E, H): dE, dH = self.computeJumps(E, H) - - if self.fluxType == "Upwind": - flux_H = 1/self.Y_imp_sum*(self.nx*self.Y_imp_p*dE-dH) - elif self.fluxType == "Centered": - flux_H = 1/self.Y_imp_sum*(self.nx*self.Y_imp_p*dE) - else: - raise ValueError("Invalid fluxType label") + flux_H = 1/self.Y_imp_sum*(self.nx*self.Y_imp_p*dE-self.alpha*dH) return flux_H - def computeFlux(self, E, H): - dE, dH = self.computeJumps(E, H) - - if self.fluxType == "Upwind": - flux_E = 1/self.Z_imp_sum*(self.nx*self.Z_imp_p*dH-dE) - flux_H = 1/self.Y_imp_sum*(self.nx*self.Y_imp_p*dE-dH) - elif self.fluxType == "Centered": - flux_E = 1/self.Z_imp_sum*(self.nx*self.Z_imp_p*dH) - flux_H = 1/self.Y_imp_sum*(self.nx*self.Y_imp_p*dE) - else: - raise ValueError("Invalid fluxType label") - return flux_E, flux_H - def computeJumps(self, E, H): Ebc, Hbc = self.fieldsOnBoundaryConditions(E, H) dE = E.transpose().take(self.vmap_m) - E.transpose().take(self.vmap_p) diff --git a/maxwell/integrators/LSERK4.py b/maxwell/integrators/LSERK4.py index 3fd59f5..61c6e60 100644 --- a/maxwell/integrators/LSERK4.py +++ b/maxwell/integrators/LSERK4.py @@ -23,6 +23,7 @@ def __init__(self, sp: SpatialDiscretization, fields): for l, f in fields.items(): self.fieldsRes[l] = np.zeros(f.shape) + def step(self, fields, dt): for s in range(0, self.N_STAGES): fieldsRHS = self.sp.computeRHS(fields) diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 16df185..b5fe3e0 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -57,34 +57,34 @@ def test_pec(): def test_pec_centered(): sp = DG1D( - n_order=5, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), + n_order=3, + mesh=Mesh1D(0, 1.0, 15, boundary_label="NULL"), fluxType="Centered" ) - driver = MaxwellDriver(sp) + driver = MaxwellDriver(sp, CFL=1.0) final_time = 1.999 - s0 = 0.25 - initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) + s0 = 0.125 + initialFieldE = np.exp(-(sp.x-0.5)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] finalFieldE = driver['E'] driver.run_until(final_time) - R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), - -finalFieldE.reshape(1, finalFieldE.size)) - assert R[0, 1] > 0.9999 + # R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), + # -finalFieldE.reshape(1, finalFieldE.size)) + # assert R[0, 1] > 0.9999 - # driver['E'][:] = initialFieldE[:] - # for _ in range(172): - # driver.step() - # plt.plot(sp.x, driver['E'],'b') - # plt.plot(sp.x, driver['H'],'r') - # plt.ylim(-1, 1) - # plt.grid(which='both') - # plt.pause(0.01) - # plt.cla() + driver['E'][:] = initialFieldE[:] + for _ in range(172): + driver.step() + plt.plot(sp.x, driver['E'],'b') + plt.plot(sp.x, driver['H'],'r') + plt.ylim(-1, 1) + plt.grid(which='both') + plt.pause(0.01) + plt.cla() def test_pec_centered_lserk74(): @@ -410,12 +410,12 @@ def test_buildCausallyConnectedOperators_2(): def test_local_causal_operator(): sp = DG1D( n_order=1, - mesh=Mesh1D(-1.0, 1.0, 11, boundary_label="Periodic"), + mesh=Mesh1D(-1.0, 1.0, 15, boundary_label="Periodic"), fluxType="Centered" ) dr = MaxwellDriver(sp) G = sp.reorder_by_elements(dr.buildDrivedEvolutionOperator()) - L = sp.reorder + L = dr.buildLocalDrivedEvolutionOperator() @@ -425,9 +425,9 @@ def test_energy_evolution_centered(): dissipate because of the LSERK4 time integration. ''' sp = DG1D( - n_order=1, + n_order=2, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp) @@ -442,26 +442,57 @@ def test_energy_evolution_centered(): driver.step() totalEnergy = energyE + energyH - # Expect some dissipation due to LSERK4. - # assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 + assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 # plt.plot(energyE, '.-b') # plt.plot(energyH, '.-r') # plt.plot(totalEnergy, '.-g') # plt.show() + + +def test_energy_evolution_upwind(): + ''' + Checks energy evolution. With Centered flux, energy should only + dissipate because of the LSERK4 time integration. + ''' + sp = DG1D( + n_order=3, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), + fluxPenalty=1.0 + ) + + driver = MaxwellDriver(sp) + driver['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) + driver['H'][:] = np.exp(-sp.x**2/(2*0.25**2)) + + Nsteps = 171 + energyE = np.zeros(Nsteps) + energyH = np.zeros(Nsteps) + for n in range(Nsteps): + energyE[n] = sp.getEnergy(driver['E']) + energyH[n] = sp.getEnergy(driver['H']) + driver.step() + + totalEnergy = energyE + energyH + assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 2e-4 + + # plt.plot(energyE, '.-b') + # plt.plot(energyH, '.-r') + # plt.plot(totalEnergy, '.-g') + plt.show() + def test_energy_evolution_with_operators(): sp = DG1D( - n_order=1, + n_order=3, mesh=Mesh1D(0, 1.0, 10, boundary_label="Periodic"), - fluxType="Upwind" - + fluxPenalty=1.0 ) dr = MaxwellDriver(sp) dr['E'][:] = np.exp(-sp.x**2/(2*0.25**2)) dr['H'][:] = np.exp(-sp.x**2/(2*0.25**2)) - + q0 = sp.fieldsAsStateVector(dr.fields) G = dr.buildDrivedEvolutionOperator() Mg = sp.buildGlobalMassMatrix() @@ -1014,24 +1045,26 @@ def test_buildDrivedEvolutionOperator(): sp = DG1D( n_order=2, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) driver = MaxwellDriver(sp) - G = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) + G = driver.buildDrivedEvolutionOperator() s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] q0 = sp.fieldsAsStateVector(driver.fields) - - driver.step() - qExpected = sp.fieldsAsStateVector(driver.fields) - - q = G.dot(q0) - - assert np.allclose(qExpected, q) + + q = np.zeros_like(q0) + q[:] = q0[:] + for _ in range(50): + driver.step() + qExpected = sp.fieldsAsStateVector(driver.fields) + q = G.dot(q) + assert np.allclose(qExpected, q) + def test_buildDrivedEvolutionOperator_LF2V(): From 6e04fb9baf95f1936ce01b46b4b3c8fbac50d17e Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Tue, 7 May 2024 22:08:18 +0200 Subject: [PATCH 10/25] fluxType changed to fluxPenalty --- maxwell/dg/dg1d.py | 3 -- test/dg/test_maxwell1d.py | 6 +-- test/test_dgtd_1d.py | 77 +++++++++++++++++---------------------- 3 files changed, 36 insertions(+), 50 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index e5f77fd..02ee516 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -88,9 +88,6 @@ def fieldsOnBoundaryConditions(self, E, H): elif label == "PMC": Hbc = - H.transpose().take(self.vmap_b) Ebc = E.transpose().take(self.vmap_b) - elif label == "NULL": - Hbc = - H.transpose().take(self.vmap_b) - Ebc = - E.transpose().take(self.vmap_b) elif label == "SMA": Hbc = H.transpose().take(self.vmap_b) * 0.0 Ebc = E.transpose().take(self.vmap_b) * 0.0 diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index 564802d..1d16d0c 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -37,7 +37,7 @@ def test_energy_with_operators(): def test_buildEvolutionOperator_PEC(): m = Mesh1D(0, 1, 5, boundary_label='PEC') - sp = DG1D(1, m, "Centered") + sp = DG1D(1, m, 0.0) A = sp.buildEvolutionOperator() A = sp.reorder_by_elements(A) M = sp.buildGlobalMassMatrix() @@ -50,7 +50,7 @@ def test_buildEvolutionOperator_PEC(): def test_buildEvolutionOperator_Periodic(): m = Mesh1D(0, 1, 5, boundary_label='Periodic') - sp = DG1D(1, m, "Centered") + sp = DG1D(1, m, 0.0) A = sp.buildEvolutionOperator() M = sp.buildGlobalMassMatrix() @@ -61,7 +61,7 @@ def test_buildEvolutionOperator_Periodic(): def test_buildEvolutionOperator_sorting(): m = Mesh1D(0, 1, 3) - sp = DG1D(2, m, "Centered") + sp = DG1D(2, m, 0.0) Np = sp.number_of_nodes_per_element() K = m.number_of_elements() diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index b5fe3e0..ce7f499 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -58,8 +58,8 @@ def test_pec(): def test_pec_centered(): sp = DG1D( n_order=3, - mesh=Mesh1D(0, 1.0, 15, boundary_label="NULL"), - fluxType="Centered" + mesh=Mesh1D(0, 1.0, 15), + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, CFL=1.0) @@ -79,19 +79,19 @@ def test_pec_centered(): driver['E'][:] = initialFieldE[:] for _ in range(172): driver.step() - plt.plot(sp.x, driver['E'],'b') - plt.plot(sp.x, driver['H'],'r') - plt.ylim(-1, 1) - plt.grid(which='both') - plt.pause(0.01) - plt.cla() + # plt.plot(sp.x, driver['E'],'b') + # plt.plot(sp.x, driver['H'],'r') + # plt.ylim(-1, 1) + # plt.grid(which='both') + # plt.pause(0.01) + # plt.cla() def test_pec_centered_lserk74(): sp = DG1D( n_order=5, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='LSERK74', CFL=1.0) @@ -127,7 +127,7 @@ def test_pec_centered_lserk134(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='LSERK134', CFL=2) @@ -158,7 +158,7 @@ def test_pec_centered_euler(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Upwind" + fluxPenalty=1.0 ) driver = MaxwellDriver(sp, timeIntegratorType='EULER', CFL=0.5659700) @@ -189,7 +189,7 @@ def test_pec_centered_lf2(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.8) @@ -220,7 +220,7 @@ def test_pec_centered_lf2v(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.8) @@ -252,7 +252,7 @@ def test_pec_centered_lf2_and_lf2v_are_equivalent(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) drLF = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.8) drVe = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.8) @@ -280,7 +280,7 @@ def test_pec_centered_ibe(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='IBE', CFL=1.5) @@ -312,7 +312,7 @@ def test_pec_centered_cn(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="PEC"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='CN', CFL=1.0) @@ -344,7 +344,7 @@ def test_periodic_centered_dirk2(): sp = DG1D( n_order=5, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) driver = MaxwellDriver(sp, timeIntegratorType='DIRK2', CFL=8) @@ -373,7 +373,7 @@ def test_buildCausallyConnectedOperators(): sp = DG1D( n_order=1, mesh=Mesh1D(-1.0, 1.0, 11, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) dr = MaxwellDriver(sp) @@ -391,7 +391,7 @@ def test_buildCausallyConnectedOperators_2(): sp = DG1D( n_order=1, mesh=Mesh1D(-1.0, 1.0, 15, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) dr = MaxwellDriver(sp) @@ -406,18 +406,7 @@ def test_buildCausallyConnectedOperators_2(): G_f = np.array([[G_reassembled[i][j] for j in new_order] for i in new_order]) assert np.all(G[:G_f.shape[0], :G_f.shape[1]] == G_f) - -def test_local_causal_operator(): - sp = DG1D( - n_order=1, - mesh=Mesh1D(-1.0, 1.0, 15, boundary_label="Periodic"), - fluxType="Centered" - ) - dr = MaxwellDriver(sp) - G = sp.reorder_by_elements(dr.buildDrivedEvolutionOperator()) - L = dr.buildLocalDrivedEvolutionOperator() - - + def test_energy_evolution_centered(): ''' @@ -479,7 +468,7 @@ def test_energy_evolution_upwind(): # plt.plot(energyE, '.-b') # plt.plot(energyH, '.-r') # plt.plot(totalEnergy, '.-g') - plt.show() + # plt.show() @@ -513,7 +502,7 @@ def test_energy_evolution_centered_lf2(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) dr = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.7) @@ -548,7 +537,7 @@ def test_energy_evolution_centered_lf2v(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) dr = MaxwellDriver(sp, timeIntegratorType='LF2V', CFL=0.7) @@ -584,7 +573,7 @@ def test_periodic_tested(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) final_time = 1.999 driver = MaxwellDriver(sp, timeIntegratorType='LSERK134') @@ -627,7 +616,7 @@ def test_periodic_LSERK_errors(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) final_time = 1.999 drLSERK54 = MaxwellDriver(sp) @@ -684,7 +673,7 @@ def test_periodic_implicit_errors(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) final_time = 1.999 @@ -726,7 +715,7 @@ def test_periodic_euler_errors(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) final_time = 1.999 @@ -780,7 +769,7 @@ def test_computational_cost(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) final_time = 1.999 drLSERK54 = MaxwellDriver(sp, CFL=20) @@ -832,7 +821,7 @@ def test_errors(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) final_time = 1.999 drLSERK54 = MaxwellDriver(sp, CFL=3) @@ -932,7 +921,7 @@ def test_max_time_step(): sp = DG1D( n_order=3, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='LSERK134') p, q = rk.loadRKM('SSP75').stability_function(mode='float') @@ -1001,7 +990,7 @@ def test_periodic_same_initial_conditions(): sp = DG1D( n_order=2, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), - fluxType="Upwind" + fluxPenalty=1.0 ) final_time = 1.999 @@ -1026,7 +1015,7 @@ def test_sma(): sp = DG1D( n_order=2, mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="SMA"), - fluxType="Upwind" + fluxPenalty=1.0 ) final_time = 3.999 @@ -1071,7 +1060,7 @@ def test_buildDrivedEvolutionOperator_LF2V(): sp = DG1D( n_order=2, mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxType="Centered" + fluxPenalty=0.0 ) driver = MaxwellDriver(sp, timeIntegratorType='LF2V') From 7d39db4b16f4a1b73be9f55b0531479840011437 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 8 May 2024 19:07:29 +0200 Subject: [PATCH 11/25] Adds ABC for centered --- maxwell/dg/dg1d.py | 47 ++++++++++++++++++++++--------------- test/test_dgtd_1d.py | 55 +++++++++++++++++++++++++++++++++----------- 2 files changed, 70 insertions(+), 32 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index 02ee516..a2c8ab7 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -7,6 +7,10 @@ class DG1D(SpatialDiscretization): def __init__(self, n_order: int, mesh: Mesh1D, fluxPenalty=1.0): + ''' + fluxPenalty is a real in [0, 1]. A value of 0.0 means no penalty + which is equivalent to a centered flux. 1.0 means full upwinding. + ''' SpatialDiscretization.__init__(self, mesh) assert n_order > 0 @@ -15,6 +19,9 @@ def __init__(self, n_order: int, mesh: Mesh1D, fluxPenalty=1.0): self.alpha = fluxPenalty self.n_faces = 2 self.n_fp = 1 + + if self.mesh.boundary_label["LEFT"] != self.mesh.boundary_label["RIGHT"]: + raise ValueError("Boundaries must be of the same type.") # Set up material parameters @@ -78,25 +85,27 @@ def get_impedance(self): return Z_imp def fieldsOnBoundaryConditions(self, E, H): - bcType = self.mesh.boundary_label - - for bdr, label in self.mesh.boundary_label.items(): - if bdr == "LEFT" or bdr == "RIGHT": - if label == "PEC": - Ebc = - E.transpose().take(self.vmap_b) - Hbc = H.transpose().take(self.vmap_b) - elif label == "PMC": - Hbc = - H.transpose().take(self.vmap_b) - Ebc = E.transpose().take(self.vmap_b) - elif label == "SMA": - Hbc = H.transpose().take(self.vmap_b) * 0.0 - Ebc = E.transpose().take(self.vmap_b) * 0.0 - elif label == "Periodic": - Ebc = E.transpose().take(self.vmap_b[::-1]) - Hbc = H.transpose().take(self.vmap_b[::-1]) - else: - raise ValueError("Invalid boundary label.") - return Ebc, Hbc + label = self.mesh.boundary_label["LEFT"] + Eb = E.transpose().take(self.vmap_b) + Hb = H.transpose().take(self.vmap_b) + if label == "PEC": + Ebc = - Eb + Hbc = Hb + elif label == "PMC": + Hbc = - Hb + Ebc = Eb + elif label == "ABC": + Ebc = (Eb + self.nx.take(self.map_b) * Hb)*0.5 + Hbc = (Hb + self.nx.take(self.map_b) * Eb)*0.5 + elif label == "SMA": + Hbc = Hb * 0.0 + Ebc = Eb * 0.0 + elif label == "Periodic": + Ebc = E.transpose().take(self.vmap_b[::-1]) + Hbc = H.transpose().take(self.vmap_b[::-1]) + else: + raise ValueError("Invalid boundary label.") + return Ebc, Hbc def computeFluxE(self, E, H): dE, dH = self.computeJumps(E, H) diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index ce7f499..70a0409 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -58,7 +58,7 @@ def test_pec(): def test_pec_centered(): sp = DG1D( n_order=3, - mesh=Mesh1D(0, 1.0, 15), + mesh=Mesh1D(0, 1.0, 15, boundary_label="PEC"), fluxPenalty=0.0 ) driver = MaxwellDriver(sp, CFL=1.0) @@ -72,19 +72,19 @@ def test_pec_centered(): driver.run_until(final_time) - # R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), - # -finalFieldE.reshape(1, finalFieldE.size)) - # assert R[0, 1] > 0.9999 + R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), + -finalFieldE.reshape(1, finalFieldE.size)) + assert R[0, 1] > 0.9999 - driver['E'][:] = initialFieldE[:] - for _ in range(172): - driver.step() - # plt.plot(sp.x, driver['E'],'b') - # plt.plot(sp.x, driver['H'],'r') - # plt.ylim(-1, 1) - # plt.grid(which='both') - # plt.pause(0.01) - # plt.cla() + # driver['E'][:] = initialFieldE[:] + # for _ in range(172): + # driver.step() + # plt.plot(sp.x, driver['E'],'b') + # plt.plot(sp.x, driver['H'],'r') + # plt.ylim(-1, 1) + # plt.grid(which='both') + # plt.pause(0.05) + # plt.cla() def test_pec_centered_lserk74(): @@ -1029,6 +1029,35 @@ def test_sma(): assert np.allclose(0.0, finalFieldE, atol=1e-6) +def test_abc_centered(): + sp = DG1D( + n_order=3, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="ABC"), + fluxPenalty=0.0 + ) + + final_time = 3.999 + driver = MaxwellDriver(sp) + initialFieldE = np.exp(-(sp.x)**2/(2*0.25**2)) + + driver['E'][:] = initialFieldE[:] + finalFieldE = driver['E'] + + driver.run_until(final_time) + + assert np.allclose(0.0, finalFieldE, atol=1e-4) + + + # driver['E'][:] = initialFieldE[:] + # for _ in range(172): + # driver.step() + # plt.plot(sp.x, driver['E'],'b') + # plt.plot(sp.x, driver['H'],'r') + # plt.ylim(-1, 1) + # plt.grid(which='both') + # plt.pause(0.05) + # plt.cla() + def test_buildDrivedEvolutionOperator(): sp = DG1D( From 09360a214e84b5228ce43a305c8f273689366c6c Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Thu, 9 May 2024 10:09:39 +0200 Subject: [PATCH 12/25] Adds stateVectorAsFields function --- maxwell/spatialDiscretization.py | 5 +++++ test/dg/test_maxwell1d.py | 22 ++++++++++++++++++++-- test/test_dgtd_1d.py | 6 +++--- 3 files changed, 28 insertions(+), 5 deletions(-) diff --git a/maxwell/spatialDiscretization.py b/maxwell/spatialDiscretization.py index 1b16227..6f9c685 100644 --- a/maxwell/spatialDiscretization.py +++ b/maxwell/spatialDiscretization.py @@ -21,6 +21,11 @@ def fieldsAsStateVector(self, fields): for f in fields.values(): q = np.append(q, f.reshape(f.size, 1, order='F')) return q + + def stateVectorAsFields(self, q): + fields = self.buildFields() + self.copyVectorToFields(q, fields) + return fields def buildStateVector(self): fields = self.buildFields() diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index 1d16d0c..a5403dd 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -81,8 +81,26 @@ def test_build_global_mass_matrix(): sp = DG1D(2, Mesh1D(0, 1, 3)) M = sp.buildGlobalMassMatrix() + N = 2 * sp.mesh.number_of_elements() * sp.number_of_nodes_per_element() + assert M.shape == (N, N) + # plt.spy(M) # plt.show() - N = 2 * sp.mesh.number_of_elements() * sp.number_of_nodes_per_element() - assert M.shape == (N, N) +def test_stateVectorAsFields(): + sp = DG1D(2, Mesh1D(0, 1, 3)) + + fields = sp.buildFields() + for l, f in fields.items(): + fields[l] = np.random.rand(*f.shape) + + q_from_fields = sp.fieldsAsStateVector(fields) + fields_from_q = sp.stateVectorAsFields(q_from_fields) + + for label, f in fields_from_q.items(): + assert np.all(fields[label] == f) + + + + + diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 70a0409..242ce8c 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -58,14 +58,14 @@ def test_pec(): def test_pec_centered(): sp = DG1D( n_order=3, - mesh=Mesh1D(0, 1.0, 15, boundary_label="PEC"), + mesh=Mesh1D(-1.0, 1.0, 15, boundary_label="PEC"), fluxPenalty=0.0 ) driver = MaxwellDriver(sp, CFL=1.0) final_time = 1.999 - s0 = 0.125 - initialFieldE = np.exp(-(sp.x-0.5)**2/(2*s0**2)) + s0 = 0.25 + initialFieldE = np.exp(-sp.x**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] finalFieldE = driver['E'] From b20ad1f5df2ab8f9cfaa37ac03a98043e91ac899 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Thu, 9 May 2024 11:05:13 +0200 Subject: [PATCH 13/25] Minor in energy --- maxwell/dg/dg1d.py | 2 +- test/dg/test_maxwell1d.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index a2c8ab7..e1ba842 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -239,7 +239,7 @@ def getEnergy(self, field): assert field.shape == (Np, K) energy = 0.0 for k in range(K): - energy += np.inner( + energy += 0.5*np.inner( field[:, k].dot(self.mass), field[:, k]*self.jacobian[:, k] ) diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index a5403dd..591f9bb 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -12,10 +12,10 @@ def test_get_energy_N1(): fields['E'].fill(0.0) fields['E'][0, 0] = 1.0 - assert np.isclose(sp.getEnergy(fields['E']), 0.1*1.0/3.0, rtol=1e-9) + assert np.isclose(sp.getEnergy(fields['E']), 0.1*1.0/3.0/2, rtol=1e-9) fields['E'].fill(1.0) - assert np.isclose(sp.getEnergy(fields['E']), 1.0, rtol=1e-9) + assert np.isclose(sp.getEnergy(fields['E']), 1.0/2, rtol=1e-9) def test_energy_with_operators(): @@ -30,7 +30,7 @@ def test_energy_with_operators(): q = sp.fieldsAsStateVector(fields) Mg = sp.buildGlobalMassMatrix() - energy = q.T.dot(Mg).dot(q) + energy = 0.5*q.T.dot(Mg).dot(q) assert np.isclose(expectedEnergy, energy) From a481ae0b4edb0d4a6e057e8ae04723d006f5ddbe Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Mon, 13 May 2024 19:40:04 +0200 Subject: [PATCH 14/25] Adds correct energy operator for FDTD. --- maxwell/dg/dg1d.py | 8 +++- maxwell/driver.py | 6 +-- maxwell/fd/fd1d.py | 71 ++++++++++++++++++++++++++++++--- test/test_dgtd_1d.py | 41 ++++++++++++++++++- test/tests_fdtd/test_fdtd_1d.py | 32 ++++++++++++++- 5 files changed, 146 insertions(+), 12 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index e1ba842..c9e9c88 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -94,6 +94,12 @@ def fieldsOnBoundaryConditions(self, E, H): elif label == "PMC": Hbc = - Hb Ebc = Eb + elif label == "Null": + Hbc = Hb + Ebc = Eb + elif label == "Double": + Hbc = -(Eb + Hb)/2 + Ebc = -(Eb - Hb)/2 elif label == "ABC": Ebc = (Eb + self.nx.take(self.map_b) * Hb)*0.5 Hbc = (Hb + self.nx.take(self.map_b) * Eb)*0.5 @@ -231,7 +237,7 @@ def buildGlobalMassMatrix(self): def getEnergy(self, field): ''' Gets energy stored in field by computing - field^T * MassMatrix * field * Jacobian. + (1/2) * field^T * MassMatrix * field * Jacobian. for each element and then the sum. ''' Np = self.number_of_nodes_per_element() diff --git a/maxwell/driver.py b/maxwell/driver.py index 96985fd..3495855 100644 --- a/maxwell/driver.py +++ b/maxwell/driver.py @@ -80,7 +80,7 @@ def run_until(self, final_time): def __getitem__(self, key): return self.fields[key] - def buildDrivedEvolutionOperator(self, reduceToEsentialDoF=True): + def buildDrivedEvolutionOperator(self, reduceToEssentialDoF=True): N = self.sp.number_of_unknowns() A = np.zeros((N,N)) @@ -95,7 +95,7 @@ def buildDrivedEvolutionOperator(self, reduceToEsentialDoF=True): self.fields = oldFields - if reduceToEsentialDoF and self.sp.isStaggered(): + if reduceToEssentialDoF and self.sp.isStaggered(): A = self.sp.reduceToEssentialDoF(A) return A @@ -115,7 +115,7 @@ def buildCausallyConnectedOperators(self, element=0, neighbors=-1): local_indices, neigh_indices = self.sp.buildLocalAndNeighborIndices(element, neighs) G = self.sp.reorder_by_elements(self.buildDrivedEvolutionOperator()) - Mg = self.sp.buildGlobalMassMatrix() + Mg = self.sp.reorder_by_elements(self.sp.buildGlobalMassMatrix()) A = G[local_indices][:,local_indices] B = G[local_indices][:,neigh_indices] C = G[neigh_indices][:,local_indices] diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index 94c4bcd..0285f73 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -6,6 +6,8 @@ from ..spatialDiscretization import * from ..dg.mesh1d import Mesh1D +import copy + class FD1D(SpatialDiscretization): def __init__(self, mesh: Mesh1D): @@ -161,6 +163,27 @@ def isStaggered(self): def number_of_nodes_per_element(self): return 1 + def number_of_unknowns(self, field='all', reduceToEssentialDoF=False): + if field == 'all': + return self.number_of_unknowns('E', reduceToEssentialDoF) \ + + self.number_of_unknowns('H', reduceToEssentialDoF) + elif field == 'E': + if reduceToEssentialDoF: + if self.mesh.boundary_label['LEFT'] == 'Periodic' and \ + self.mesh.boundary_label['RIGHT'] == 'Periodic': + return len(self.x) - 1 + elif self.mesh.boundary_label['LEFT'] == 'PEC' and \ + self.mesh.boundary_label['LEFT'] == 'PEC': + return len(self.x) - 2 + else: + raise ValueError('Invalid boundary labels for reduction.') + else: + return len(self.x) + elif field == 'H': + return len(self.xH) + else: + raise ValueError('Invalid field label.') + def setFieldWithIndex(self, fields, i, val): NE = fields['E'].size if i < NE: @@ -180,12 +203,12 @@ def reduceToEssentialDoF(self, A): and self.mesh.boundary_label['RIGHT'] == 'PEC': A = np.delete(A, NE-1, 0) A = np.delete(A, NE-1, 1) - A = np.delete(A, 0, 0) + A = np.delete(A, 0, 0) A = np.delete(A, 0, 1) else: raise ValueError( "Periodic conditions must be ensured at both ends") - + return A def buildEvolutionOperator(self, reduceToEssentialDoF=True): @@ -197,11 +220,47 @@ def buildEvolutionOperator(self, reduceToEssentialDoF=True): fieldsRHS = self.computeRHS(fields) q0 = np.concatenate([fieldsRHS['E'], fieldsRHS['H']]) A[:, i] = q0[:] - + if reduceToEssentialDoF: - A = self.reduceToEssentialDoF(A) + A = self.reduceToEssentialDoF(A) return A + def getEnergy(self, field): + h = self.x[1] - self.x[0] + assert np.allclose(h, self.x[1:] - self.x[:-1]) + + M = np.eye(len(field)) * h + return 0.5 * field.T.dot(M).dot(field) + + def getTotalEnergy(self, G, fields): + dt = self.dt + + Gdt = G*dt + + N = self.number_of_unknowns( reduceToEssentialDoF=True) + NE = self.number_of_unknowns('E', reduceToEssentialDoF=True) + NH = self.number_of_unknowns('H', reduceToEssentialDoF=True) + + S_E = np.zeros((N, N)) + S_E[:NE, :NE] = np.eye(NE) + S_H = np.zeros((N, N)) + S_H[NE:, NE:] = np.eye(NH) + + h = self.x[1] - self.x[0] + assert np.allclose(h, self.x[1:] - self.x[:-1]) + + M = np.eye(N)*h + V = 0.5*(M - 0.5*S_E.T.dot(Gdt.T).dot(S_H) - 0.5*S_H.T.dot(Gdt).dot(S_E)) + + if self.mesh.boundary_label['LEFT'] != 'Periodic': + raise ValueError("Only implemented for periodic.") + + f = copy.deepcopy(fields) + f['E'] = np.zeros(len(fields['E'])-1) + f['E'][:] = fields['E'][:-1] + q = self.fieldsAsStateVector(f) + return q.T.dot(V).dot(q) + def reorder_by_elements(self, A): # Assumes that the original array contains all DoF ordered as: # [ E_0, ..., E_{NE-1}, H_0, ..., H_{NH-1} ] @@ -219,13 +278,13 @@ def reorder_by_elements(self, A): "Unable to order by elements with different size fields.") N = NE + NH new_order = np.zeros(N, dtype=int) - 1 - + for i in range(N): if i < NE: new_order[2*i] = i else: new_order[2*int(i - NE)+1] = i - + if (len(A.shape) == 1): A1 = [A[i] for i in new_order] elif (len(A.shape) == 2): diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 242ce8c..08197ee 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -470,6 +470,45 @@ def test_energy_evolution_upwind(): # plt.plot(totalEnergy, '.-g') # plt.show() +def test_abnormal_energy_evolution_upwind(): + ''' + This abnormal "energy" evolution happens because the upwind fluxes + are not accounted correctly to compute the physical energy. + ''' + sp = DG1D( + n_order=3, + mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), + fluxPenalty=1.0 + ) + dr = MaxwellDriver(sp) + + P = dr.buildPowerOperator() + wP, vP = np.linalg.eig(P) + + unstableIndices = np.where(wP>0.0)[0] + assert len(unstableIndices) > 0 + + dr.fields = sp.stateVectorAsFields(np.real(vP[:,unstableIndices[0]])) + + Nsteps = 5 + energyE = np.zeros(Nsteps) + energyH = np.zeros(Nsteps) + for n in range(Nsteps): + energyE[n] = sp.getEnergy(dr['E']) + energyH[n] = sp.getEnergy(dr['H']) + # plt.plot(sp.x, dr['E'],'b') + # plt.plot(sp.x, dr['H'],'r') + # plt.show() + dr.step() + + totalEnergy = energyE + energyH + assert (totalEnergy[1]-totalEnergy[0]) > 0 + assert np.all(totalEnergy[2:]-totalEnergy[1:-1] < 0) + + # plt.plot(energyE, '.-b') + # plt.plot(energyH, '.-r') + plt.plot(totalEnergy, '.-g') + plt.show() def test_energy_evolution_with_operators(): @@ -1093,7 +1132,7 @@ def test_buildDrivedEvolutionOperator_LF2V(): ) driver = MaxwellDriver(sp, timeIntegratorType='LF2V') - G = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) + G = driver.buildDrivedEvolutionOperator(reduceToEssentialDoF=False) s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index 4a36416..8f0297d 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -27,7 +27,7 @@ def test_buildDrivedEvolutionOperator(): sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="PEC")) driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) - A = driver.buildDrivedEvolutionOperator(reduceToEsentialDoF=False) + A = driver.buildDrivedEvolutionOperator(reduceToEssentialDoF=False) s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) @@ -218,6 +218,36 @@ def test_fdtd_check_initial_conditions_GW_right(): R1 = np.corrcoef(expectedE, evolvedE) assert R1[0, 1] > 0.995 + +def test_energy_evolution(): + ''' + Eneregy evolution for LF2 needs to account for the fact that + the magnetic field is staggered in time. + This requires a special operator to compute energy. + ''' + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="Periodic")) + dr = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) + + G = dr.buildDrivedEvolutionOperator(reduceToEssentialDoF=True) + s0 = 0.25 + dr['E'][:] = np.exp(-(sp.x)**2/(2*s0**2)) + + + Nsteps = 500 + energyE = np.zeros(Nsteps) + energyH = np.zeros(Nsteps) + totalEnergy = np.zeros(Nsteps) + for n in range(Nsteps): + energyE[n] = sp.getEnergy(dr['E']) + energyH[n] = sp.getEnergy(dr['H']) + totalEnergy[n] = sp.getTotalEnergy(G, dr.fields) + dr.step() + + # plt.plot(energyE + energyH) + # plt.plot((energyE[:-1] +energyE[1:])*0.5 + energyH[:-1]) + # plt.plot(totalEnergy) + # plt.show() + assert np.isclose(totalEnergy[0],totalEnergy[-1]) def test_fdtd_periodic_lserk(): From 4fdcd09b34ef0ff48418299771f70b5625ce204e Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Fri, 17 May 2024 11:53:31 +0200 Subject: [PATCH 15/25] Found bug in driven evolution operator --- maxwell/fd/fd1d.py | 12 ++++++++---- maxwell/integrators/LF2.py | 4 ++-- test/tests_fdtd/test_fdtd_1d.py | 14 ++++++++------ 3 files changed, 18 insertions(+), 12 deletions(-) diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index 0285f73..a141ee5 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -228,9 +228,13 @@ def buildEvolutionOperator(self, reduceToEssentialDoF=True): def getEnergy(self, field): h = self.x[1] - self.x[0] assert np.allclose(h, self.x[1:] - self.x[:-1]) - - M = np.eye(len(field)) * h - return 0.5 * field.T.dot(M).dot(field) + f = copy.deepcopy(field) + if self.mesh.boundary_label['LEFT'] == 'Periodic': + f = np.zeros(len(field)-1) + f[:] = field[:-1] + + M = np.eye(len(f)) * h + return 0.5 * f.T.dot(M).dot(f) def getTotalEnergy(self, G, fields): dt = self.dt @@ -250,7 +254,7 @@ def getTotalEnergy(self, G, fields): assert np.allclose(h, self.x[1:] - self.x[:-1]) M = np.eye(N)*h - V = 0.5*(M - 0.5*S_E.T.dot(Gdt.T).dot(S_H) - 0.5*S_H.T.dot(Gdt).dot(S_E)) + V = 0.5*(M + 0.5*S_E.T.dot(Gdt.T).dot(S_H) + 0.5*S_H.T.dot(Gdt).dot(S_E)) if self.mesh.boundary_label['LEFT'] != 'Periodic': raise ValueError("Only implemented for periodic.") diff --git a/maxwell/integrators/LF2.py b/maxwell/integrators/LF2.py index 1482696..4c08165 100644 --- a/maxwell/integrators/LF2.py +++ b/maxwell/integrators/LF2.py @@ -13,10 +13,10 @@ def step(self, fields, dt): H = fields['H'] if self.sp.dimension() == 1: - self.time += dt/2 - E += dt * self.sp.computeRHSE(fields) self.time += dt/2 H += dt * self.sp.computeRHSH(fields) + self.time += dt/2 + E += dt * self.sp.computeRHSE(fields) elif self.sp.dimension() == 2: self.time += dt/2 rhsE = self.sp.computeRHSE(fields) diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index 8f0297d..610a8ee 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -71,7 +71,8 @@ def test_fdtd_pec(): initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] - # plot(sp, driver) + plot(sp, driver) + driver.run_until(2.0) @@ -81,7 +82,8 @@ def test_fdtd_pec(): def test_fdtd_periodic(): - sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 5, boundary_label="Periodic")) + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="Periodic")) + driver = MaxwellDriver(sp, timeIntegratorType='LF2') s0 = 0.25 @@ -243,10 +245,10 @@ def test_energy_evolution(): totalEnergy[n] = sp.getTotalEnergy(G, dr.fields) dr.step() - # plt.plot(energyE + energyH) - # plt.plot((energyE[:-1] +energyE[1:])*0.5 + energyH[:-1]) - # plt.plot(totalEnergy) - # plt.show() + plt.plot(energyE + energyH) + plt.plot((energyE[:-1] + energyE[1:])*0.5 + energyH[:-1]) + plt.plot(totalEnergy) + plt.show() assert np.isclose(totalEnergy[0],totalEnergy[-1]) From 1a9f9358016f62db9663b0078468329cbd4aabc9 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 22 May 2024 13:36:52 +0200 Subject: [PATCH 16/25] Fixes error in FDTD energy calculation for periodic. --- maxwell/fd/fd1d.py | 7 +++---- test/test_dgtd_1d.py | 25 +++++++++++++------------ test/tests_fdtd/test_fdtd_1d.py | 21 ++++++++++++--------- 3 files changed, 28 insertions(+), 25 deletions(-) diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index a141ee5..861d025 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -225,16 +225,15 @@ def buildEvolutionOperator(self, reduceToEssentialDoF=True): A = self.reduceToEssentialDoF(A) return A - def getEnergy(self, field): + def getEnergy(self, field, removeLast=False): h = self.x[1] - self.x[0] assert np.allclose(h, self.x[1:] - self.x[:-1]) f = copy.deepcopy(field) - if self.mesh.boundary_label['LEFT'] == 'Periodic': + if removeLast: f = np.zeros(len(field)-1) f[:] = field[:-1] - M = np.eye(len(f)) * h - return 0.5 * f.T.dot(M).dot(f) + return 0.5 * h * f.T.dot(f) def getTotalEnergy(self, G, fields): dt = self.dt diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 08197ee..7917f56 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -507,8 +507,8 @@ def test_abnormal_energy_evolution_upwind(): # plt.plot(energyE, '.-b') # plt.plot(energyH, '.-r') - plt.plot(totalEnergy, '.-g') - plt.show() + # plt.plot(totalEnergy, '.-g') + # plt.show() def test_energy_evolution_with_operators(): @@ -536,11 +536,11 @@ def test_energy_evolution_with_operators(): assert np.isclose(power, expectedPower) - +@pytest.mark.skip(reason="Unsure about the result.") def test_energy_evolution_centered_lf2(): sp = DG1D( n_order=3, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), + mesh=Mesh1D(-5.0, 5.0, 50, boundary_label="Periodic"), fluxPenalty=0.0 ) @@ -561,21 +561,22 @@ def test_energy_evolution_centered_lf2(): # plt.cla() dr.step() - totalEnergy = energyE[1:] + (energyH[:-1]+ energyH[1:])*.5 + totalEnergy = energyE[:1] + (energyH[:-1]+ energyH[1:])*.5 + - assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 1e-6 + # assert np.abs(totalEnergy[-1]-totalEnergy[0]) < 1e-6 - # plt.plot(energyE) - # plt.plot(energyH) - # plt.plot(totalEnergy, '-b') - # plt.grid() - # plt.show() + plt.plot(energyE) + plt.plot(energyH) + plt.plot(totalEnergy, '-b') + plt.grid() + plt.show() def test_energy_evolution_centered_lf2v(): sp = DG1D( n_order=3, - mesh=Mesh1D(-1.0, 1.0, 10, boundary_label="Periodic"), + mesh=Mesh1D(-5.0, 5.0, 50, boundary_label="Periodic"), fluxPenalty=0.0 ) diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index 610a8ee..8cdae08 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -71,8 +71,7 @@ def test_fdtd_pec(): initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) driver['E'][:] = initialFieldE[:] - plot(sp, driver) - + # plot(sp, driver) driver.run_until(2.0) @@ -223,7 +222,7 @@ def test_fdtd_check_initial_conditions_GW_right(): def test_energy_evolution(): ''' - Eneregy evolution for LF2 needs to account for the fact that + Energy evolution for LF2 needs to account for the fact that the magnetic field is staggered in time. This requires a special operator to compute energy. ''' @@ -234,21 +233,25 @@ def test_energy_evolution(): s0 = 0.25 dr['E'][:] = np.exp(-(sp.x)**2/(2*s0**2)) + if sp.mesh.boundary_label['LEFT'] == 'Periodic': + removeLastE = True - Nsteps = 500 + Nsteps = 300 energyE = np.zeros(Nsteps) energyH = np.zeros(Nsteps) totalEnergy = np.zeros(Nsteps) for n in range(Nsteps): - energyE[n] = sp.getEnergy(dr['E']) + energyE[n] = sp.getEnergy(dr['E'], removeLast=removeLastE) energyH[n] = sp.getEnergy(dr['H']) totalEnergy[n] = sp.getTotalEnergy(G, dr.fields) dr.step() - plt.plot(energyE + energyH) - plt.plot((energyE[:-1] + energyE[1:])*0.5 + energyH[:-1]) - plt.plot(totalEnergy) - plt.show() + # plt.plot(energyE + energyH) + # plt.plot((energyH[:-1] + energyH[1:])*0.5 + energyE[:-1]) + # plt.plot(energyE) + # plt.plot(energyH) + # plt.plot(totalEnergy) + # plt.show() assert np.isclose(totalEnergy[0],totalEnergy[-1]) From 30513f405027745fa1d7438d0df62f40ac349d7e Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Thu, 23 May 2024 10:40:09 +0200 Subject: [PATCH 17/25] Fixes Mur FDTD in 1D. Was failing for CFL < 1.0 --- maxwell/fd/fd1d.py | 27 ++++++++++++++++----------- test/tests_fdtd/test_fdtd_1d.py | 12 ++++++------ 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index 861d025..c93cd7b 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -76,43 +76,48 @@ def computeRHSE(self, fields): if label == "PEC": rhsE[0] = 0.0 - if label == "PMC": + elif label == "PMC": rhsE[0] = - (1.0/self.dxH[0]) * (2 * H[0]) - if label == "Periodic": + elif label == "Periodic": rhsE[0] = - (1.0/self.dxH[0]) * (H[0] - H[-1]) rhsE[-1] = rhsE[0] - if label == "Mur": - + elif label == "Mur": rhsE[0] = E[1] + \ (self.c0 * self.dt - self.dx[0]) / \ (self.c0 * self.dt + self.dx[0]) * \ - (rhsE[1] - E[0]) + (rhsE[1]*self.dt + E[1] - E[0]) rhsE[0] -= E[0] rhsE[0] /= self.dt + + else: + raise ValueError("Invalid boundary.") if bdr == "RIGHT": if label == "PEC": rhsE[-1] = 0.0 - if label == "PMC": + elif label == "PMC": rhsE[-1] = - (1.0/self.dxH[0]) * (-2 * H[-1]) - if label == "Periodic": + elif label == "Periodic": rhsE[0] = - (1.0/self.dxH[0]) * (H[0] - H[-1]) rhsE[-1] = rhsE[0] - if label == "Mur": + elif label == "Mur": rhsE[-1] = E[-2] + \ - (self.c0 * self.dt - self.dx[0]) / \ - (self.c0 * self.dt + self.dx[0]) * \ - (rhsE[-2] - E[-1]) + (self.c0 * self.dt - self.dx[-1]) / \ + (self.c0 * self.dt + self.dx[-1]) * \ + (rhsE[-2]*self.dt + E[-2] - E[-1]) rhsE[-1] -= E[-1] rhsE[-1] /= self.dt + + else: + raise ValueError("Invalid boundary.") return rhsE diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index 8cdae08..d0508a1 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -136,7 +136,7 @@ def test_fdtd_pmc_cfl_equals_half(): def test_fdtd_mur(): sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="Mur")) - driver = MaxwellDriver(sp, timeIntegratorType='LF2') + driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) s0 = 0.25 initialFieldE = np.exp(-(sp.x)**2/(2*s0**2)) @@ -159,7 +159,7 @@ def test_fdtd_mur_right_only(): s0 = 0.25 driver['E'][:] = np.exp(-(sp.x)**2/(2*s0**2)) - driver['H'][:] = np.exp(-(sp.xH - driver.dt/2)**2/(2*s0**2)) + driver['H'][:] = np.exp(-(sp.xH + driver.dt/2)**2/(2*s0**2)) # plot(sp, driver) @@ -178,13 +178,13 @@ def test_fdtd_right_only_mur_and_pec(): t_final = 8.0 - sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label=bdrs)) - driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) - + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 50, boundary_label=bdrs)) + driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.7) + s0 = 0.25 driver['E'][:] = np.exp(-(sp.x)**2/(2*s0**2)) initialFieldE = driver['E'][:] - driver['H'][:] = np.exp(-(sp.xH - driver.dt/2)**2/(2*s0**2)) + driver['H'][:] = -np.exp(-(sp.xH - driver.dt/2)**2/(2*s0**2)) # plot(sp, driver) From 084c8e2f0fcfe8413411b160fd7ab48faa74ce71 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Thu, 23 May 2024 13:48:08 +0200 Subject: [PATCH 18/25] Mur constant mode --- maxwell/fd/fd1d.py | 2 ++ test/tests_fdtd/test_fdtd_1d.py | 16 ++++++++++++++++ 2 files changed, 18 insertions(+) diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index c93cd7b..03aa5eb 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -75,6 +75,7 @@ def computeRHSE(self, fields): if bdr == "LEFT": if label == "PEC": rhsE[0] = 0.0 + # rhsE[0] = (0.0 - E[0])/self.dt elif label == "PMC": rhsE[0] = - (1.0/self.dxH[0]) * (2 * H[0]) @@ -98,6 +99,7 @@ def computeRHSE(self, fields): if bdr == "RIGHT": if label == "PEC": rhsE[-1] = 0.0 + #rhsE[-1] = (0.0 - E[-1])/self.dt elif label == "PMC": rhsE[-1] = - (1.0/self.dxH[0]) * (-2 * H[-1]) diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index d0508a1..810370d 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -148,6 +148,22 @@ def test_fdtd_mur(): finalFieldE = driver['E'][:] assert np.allclose(finalFieldE, 0.0, atol=1e-3) + +def test_fdtd_mur_constant_mode(): + sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label='Mur')) + + driver = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) + + s0 = 0.25 + driver['E'][:] = 0.1 + driver['H'][:] = 0.1 + + # plot(sp, driver) + + driver.run_until(1.0) + + finalFieldE = driver['E'][:] + assert np.allclose(finalFieldE, 0.1, atol=1e-12) def test_fdtd_mur_right_only(): From 1691a22e17a2401d42285e48a1c567e768ed0a82 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Fri, 24 May 2024 19:38:56 +0200 Subject: [PATCH 19/25] extends number of unknowns for dgtd --- maxwell/dg/dg1d.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index c9e9c88..1e0dbd3 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -221,6 +221,14 @@ def reorder_by_elements(self, A): A1 = [[A[i][j] for j in new_order] for i in new_order] return np.array(A1) + def number_of_unknowns(self, field='all'): + if field == 'all': + return self.number_of_unknowns('E') \ + + self.number_of_unknowns('H') + else: + f = self.buildFields() + return f[field].size + def buildGlobalMassMatrix(self): Np = self.number_of_nodes_per_element() K = self.mesh.number_of_elements() From 8936e56e96079f99a96a4e2bd18248ab2c18a2b0 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Mon, 27 May 2024 18:11:32 +0200 Subject: [PATCH 20/25] Flux and Stiffness operators --- maxwell/dg/dg1d.py | 44 ++++++++++++++++++++++++++++++++++ test/dg/test_maxwell1d.py | 12 ++++++++++ test/test_dgtd_1d.py | 50 +++++++++++++++++++++++++++++++++++---- 3 files changed, 101 insertions(+), 5 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index 1e0dbd3..d11aaf7 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -229,6 +229,50 @@ def number_of_unknowns(self, field='all'): f = self.buildFields() return f[field].size + def buildStiffnessMatrix(self): + Np = self.number_of_nodes_per_element() + K = self.mesh.number_of_elements() + N = self.number_of_unknowns() + A = np.zeros((N, N)) + for i in range(N): + fields = self.buildFields() + E = fields['E'] + H = fields['H'] + self.setFieldWithIndex(fields, i, 1.0) + rhs_drE = np.matmul(self.diff_matrix, E) + rhs_drH = np.matmul(self.diff_matrix, H) + rhsE = 1/self.epsilon * np.multiply(-1*self.rx, rhs_drH) + rhsH = 1/self.mu * np.multiply(-1*self.rx, rhs_drE) + q0 = np.vstack([ + rhsE.reshape(Np*K, 1, order='F'), + rhsH.reshape(Np*K, 1, order='F') + ]) + A[:, i] = q0[:, 0] + return A + + + + def buildFluxMatrix(self): + Np = self.number_of_nodes_per_element() + K = self.mesh.number_of_elements() + N = self.number_of_unknowns() + A = np.zeros((N, N)) + for i in range(N): + fields = self.buildFields() + E = fields['E'] + H = fields['H'] + self.setFieldWithIndex(fields, i, 1.0) + flux_E = self.computeFluxE(E, H) + flux_H = self.computeFluxH(E, H) + rhsE = 1/self.epsilon * np.matmul(self.lift, self.f_scale * flux_E) + rhsH = 1/self.mu * (np.matmul(self.lift, self.f_scale * flux_H)) + q0 = np.vstack([ + rhsE.reshape(Np*K, 1, order='F'), + rhsH.reshape(Np*K, 1, order='F') + ]) + A[:, i] = q0[:, 0] + return A + def buildGlobalMassMatrix(self): Np = self.number_of_nodes_per_element() K = self.mesh.number_of_elements() diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index 591f9bb..ee137ca 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -57,7 +57,19 @@ def test_buildEvolutionOperator_Periodic(): assert np.allclose(A.T.dot(M) + (M).dot(A), 0.0) assert A.shape == (20, 20) assert np.allclose(np.real(np.linalg.eig(A)[0]), 0) + + +def test_stiffness_and_flux_operators(): + m = Mesh1D(0, 1, 5, boundary_label='Periodic') + sp = DG1D(1, m, 0.0) + + A = sp.buildEvolutionOperator() + S = sp.buildStiffnessMatrix() + F = sp.buildFluxMatrix() + + assert np.allclose(A, S+F) + def test_buildEvolutionOperator_sorting(): m = Mesh1D(0, 1, 3) diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index 7917f56..c978e0b 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -340,6 +340,46 @@ def test_pec_centered_cn(): # plt.cla() + +def test_periodic_centered(): + sp = DG1D( + n_order=3, + mesh=Mesh1D(0, 1.0, 30, boundary_label="Periodic"), + fluxPenalty=0.0 + ) + dr = MaxwellDriver(sp, timeIntegratorType='LSERK4', CFL=1.365) + + # final_time = 1.999 + s0 = 0.125 + initialFieldE = np.exp(-(sp.x-0.5)**2/(2*s0**2)) + + # dr['E'][:] = initialFieldE[:] + # finalFieldE = dr['E'] + + # dr.run_until(final_time) + + # R = np.corrcoef(initialFieldE.reshape(1, initialFieldE.size), + # -finalFieldE.reshape(1, finalFieldE.size)) + # assert R[0, 1] > 0.9999 + + + dr['E'][:] = initialFieldE[:] + + initialEnergy = sp.getEnergy(dr['E']) + sp.getEnergy(dr['H']) + for _ in range(500): + dr.step() + # plt.plot(sp.x, dr['E'],'b') + # plt.plot(sp.x, dr['H'],'r') + # plt.ylim(-1, 1) + # plt.grid(which='both') + # plt.pause(0.05) + # plt.cla() + finalEnergy = sp.getEnergy(dr['E']) + sp.getEnergy(dr['H']) + maxEVG = np.max(np.abs(np.linalg.eig(dr.buildDrivedEvolutionOperator())[0])) + + assert maxEVG < 1.0 or np.isclose(maxEVG, 1.0) + assert finalEnergy <= initialEnergy + def test_periodic_centered_dirk2(): sp = DG1D( n_order=5, @@ -612,11 +652,11 @@ def test_energy_evolution_centered_lf2v(): def test_periodic_tested(): sp = DG1D( n_order=3, - mesh=Mesh1D(-1.0, 1.0, 20, boundary_label="Periodic"), - fluxPenalty=1.0 + mesh=Mesh1D(-1.0, 1.0, 30, boundary_label="Periodic"), + fluxPenalty=0.0 ) final_time = 1.999 - driver = MaxwellDriver(sp, timeIntegratorType='LSERK134') + driver = MaxwellDriver(sp, timeIntegratorType='LSERK54') initialField = np.sin(2*np.pi*sp.x) driver['E'][:] = initialField[:] @@ -648,7 +688,7 @@ def real_function(x, t): driver.step() t += driver.dt - assert (np.sqrt(error).max() < 1e-02, True) + assert np.sqrt(error).max() < 1e-08 @pytest.mark.skip(reason="Nothing is being tested.") @@ -953,7 +993,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.sqrt(error_abs).max() < 1e-02 @pytest.mark.skip(reason="Nothing is being tested.") From 6c492e5cf12be01d36c5b2ad100260d14bc509c5 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Fri, 31 May 2024 13:55:29 +0200 Subject: [PATCH 21/25] [WIP] SP connected operators --- maxwell/dg/dg1d.py | 18 ++++++++++++++ test/dg/test_dg1d.py | 1 + test/dg/test_maxwell1d.py | 51 +++++++++++++++++++++++---------------- test/test_dgtd_1d.py | 2 +- 4 files changed, 50 insertions(+), 22 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index d11aaf7..500b03a 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -291,6 +291,7 @@ def getEnergy(self, field): Gets energy stored in field by computing (1/2) * field^T * MassMatrix * field * Jacobian. for each element and then the sum. + WRONG IF EVOLVED BY A STAGGERED SCHEME. ''' Np = self.number_of_nodes_per_element() K = self.mesh.number_of_elements() @@ -303,3 +304,20 @@ def getEnergy(self, field): ) return energy + + def buildConnectedOperators(self, element=0, neighbors=1): + + neighs = neighbors + + local_indices, neigh_indices = self.buildLocalAndNeighborIndices(element, neighs) + + G = self.reorder_by_elements(self.buildEvolutionOperator()) + Mg = self.reorder_by_elements(self.buildGlobalMassMatrix()) + A = G[local_indices][:,local_indices] + B = G[local_indices][:,neigh_indices] + C = G[neigh_indices][:,local_indices] + D = G[neigh_indices][:,neigh_indices] + Mk = Mg[local_indices][:,local_indices] + Mn = Mg[neigh_indices][:,neigh_indices] + + return A, B, C, D, Mk, Mn \ No newline at end of file diff --git a/test/dg/test_dg1d.py b/test/dg/test_dg1d.py index 6b7aaac..e3ca229 100644 --- a/test/dg/test_dg1d.py +++ b/test/dg/test_dg1d.py @@ -254,3 +254,4 @@ def test_mass_matrix(): [-1., 2., 4.]]) / 15 assert np.allclose(M, expected) + diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index ee137ca..852b7da 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -21,19 +21,19 @@ def test_get_energy_N1(): def test_energy_with_operators(): m = Mesh1D(0, 1, 10) sp = DG1D(1, m) - + fields = sp.buildFields() fields['E'].fill(1.0) fields['H'].fill(2.0) - + expectedEnergy = sp.getEnergy(fields['E']) + sp.getEnergy(fields['H']) - + q = sp.fieldsAsStateVector(fields) Mg = sp.buildGlobalMassMatrix() energy = 0.5*q.T.dot(Mg).dot(q) - + assert np.isclose(expectedEnergy, energy) - + def test_buildEvolutionOperator_PEC(): m = Mesh1D(0, 1, 5, boundary_label='PEC') @@ -42,7 +42,6 @@ def test_buildEvolutionOperator_PEC(): A = sp.reorder_by_elements(A) M = sp.buildGlobalMassMatrix() - assert np.allclose(A.T.dot(M) + (M).dot(A), 0.0) assert A.shape == (20, 20) assert np.allclose(np.real(np.linalg.eig(A)[0]), 0) @@ -57,19 +56,18 @@ def test_buildEvolutionOperator_Periodic(): assert np.allclose(A.T.dot(M) + (M).dot(A), 0.0) assert A.shape == (20, 20) assert np.allclose(np.real(np.linalg.eig(A)[0]), 0) - - + + def test_stiffness_and_flux_operators(): m = Mesh1D(0, 1, 5, boundary_label='Periodic') sp = DG1D(1, m, 0.0) - + A = sp.buildEvolutionOperator() S = sp.buildStiffnessMatrix() F = sp.buildFluxMatrix() - + assert np.allclose(A, S+F) - def test_buildEvolutionOperator_sorting(): m = Mesh1D(0, 1, 3) @@ -89,30 +87,41 @@ def test_buildEvolutionOperator_sorting(): np.sort(np.imag(eigA_by_elem))) +def test_build_connected_operators(): + sp = DG1D(2, Mesh1D(0, 1, 3), 0.0) + + Ag = sp.reorder_by_elements(sp.buildEvolutionOperator()) + eigAg = np.sort(np.linalg.eig(Ag)[0]) + A, B, C, D, _, _ = sp.buildConnectedOperators(1, 1) + + Ag_reassembled_1 = np.concatenate([A, B], axis=1) + Ag_reassembled_2 = np.concatenate([C, D], axis=1) + Ag_reassembled = np.concatenate([Ag_reassembled_1, Ag_reassembled_2]) + eigAg_reassembled = np.sort(np.linalg.eig(Ag_reassembled)[0]) + + assert np.allclose(eigAg, eigAg_reassembled) + + def test_build_global_mass_matrix(): sp = DG1D(2, Mesh1D(0, 1, 3)) M = sp.buildGlobalMassMatrix() N = 2 * sp.mesh.number_of_elements() * sp.number_of_nodes_per_element() assert M.shape == (N, N) - + # plt.spy(M) # plt.show() + def test_stateVectorAsFields(): sp = DG1D(2, Mesh1D(0, 1, 3)) - + fields = sp.buildFields() for l, f in fields.items(): fields[l] = np.random.rand(*f.shape) - + q_from_fields = sp.fieldsAsStateVector(fields) fields_from_q = sp.stateVectorAsFields(q_from_fields) - + for label, f in fields_from_q.items(): - assert np.all(fields[label] == f) - - - - - + assert np.all(fields[label] == f) diff --git a/test/test_dgtd_1d.py b/test/test_dgtd_1d.py index c978e0b..3294590 100644 --- a/test/test_dgtd_1d.py +++ b/test/test_dgtd_1d.py @@ -648,7 +648,7 @@ def test_energy_evolution_centered_lf2v(): # plt.grid() # plt.show() - +@pytest.mark.skip(reason="Unsure about the result.") def test_periodic_tested(): sp = DG1D( n_order=3, From 76cb5a24ce499920781a26aeeca8ea721c4192cf Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Sat, 1 Jun 2024 13:36:32 +0200 Subject: [PATCH 22/25] Minor --- test/dg/test_maxwell1d.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/test/dg/test_maxwell1d.py b/test/dg/test_maxwell1d.py index 852b7da..0c4d94a 100644 --- a/test/dg/test_maxwell1d.py +++ b/test/dg/test_maxwell1d.py @@ -88,19 +88,25 @@ def test_buildEvolutionOperator_sorting(): def test_build_connected_operators(): - sp = DG1D(2, Mesh1D(0, 1, 3), 0.0) - + sp = DG1D(2, Mesh1D(0, 1, 11, "Periodic"), 0.0) + + Np = sp.number_of_nodes_per_element() + k = 5 + n = 1 + Ag = sp.reorder_by_elements(sp.buildEvolutionOperator()) - eigAg = np.sort(np.linalg.eig(Ag)[0]) - A, B, C, D, _, _ = sp.buildConnectedOperators(1, 1) - - Ag_reassembled_1 = np.concatenate([A, B], axis=1) - Ag_reassembled_2 = np.concatenate([C, D], axis=1) - Ag_reassembled = np.concatenate([Ag_reassembled_1, Ag_reassembled_2]) - eigAg_reassembled = np.sort(np.linalg.eig(Ag_reassembled)[0]) - - assert np.allclose(eigAg, eigAg_reassembled) - + A, B, C, D, _, _ = sp.buildConnectedOperators(k, n) + + assert np.allclose(A, Ag[k*2*Np:(k+1)*2*Np,k*2*Np:(k+1)*2*Np]) + assert np.allclose(B[:, :n*2*Np], Ag[k*2*Np:(k+1)*2*Np,(k-n)*2*Np:k*2*Np]) + assert np.allclose(B[:, n*2*Np:], Ag[k*2*Np:(k+1)*2*Np,(k+1)*2*Np:(k+1+n)*2*Np]) + assert np.allclose(C[:n*2*Np, :], Ag[(k-n)*2*Np:k*2*Np,k*2*Np:(k+1)*2*Np]) + assert np.allclose(C[n*2*Np:, :], Ag[(k+1)*2*Np:(k+1+n)*2*Np,k*2*Np:(k+1)*2*Np]) + assert np.allclose(D[:2*Np, :2*Np], A) + assert np.allclose(D[2*Np:, 2*Np:], A) + assert np.allclose(D[:2*Np, 2*Np:], 0) + assert np.allclose(D[2*Np:, :2*Np], 0) + def test_build_global_mass_matrix(): sp = DG1D(2, Mesh1D(0, 1, 3)) From 28b354b4e93ee6c076c2b508d9956cfee5dbe8cb Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Wed, 14 Aug 2024 13:01:29 +0200 Subject: [PATCH 23/25] Minor in fdtd total energy calculation --- maxwell/fd/fd1d.py | 19 ++++++++----------- test/tests_fdtd/test_fdtd_1d.py | 6 +++--- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index 03aa5eb..abb0e43 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -243,24 +243,21 @@ def getEnergy(self, field, removeLast=False): return 0.5 * h * f.T.dot(f) def getTotalEnergy(self, G, fields): - dt = self.dt - - Gdt = G*dt - N = self.number_of_unknowns( reduceToEssentialDoF=True) NE = self.number_of_unknowns('E', reduceToEssentialDoF=True) NH = self.number_of_unknowns('H', reduceToEssentialDoF=True) - S_E = np.zeros((N, N)) - S_E[:NE, :NE] = np.eye(NE) - S_H = np.zeros((N, N)) - S_H[NE:, NE:] = np.eye(NH) + L_E = np.zeros((N, N)) + L_E[:NE, :NE] = np.eye(NE) + L_H = np.zeros((N, N)) + L_H[NE:, NE:] = np.eye(NH) h = self.x[1] - self.x[0] assert np.allclose(h, self.x[1:] - self.x[:-1]) - M = np.eye(N)*h - V = 0.5*(M + 0.5*S_E.T.dot(Gdt.T).dot(S_H) + 0.5*S_H.T.dot(Gdt).dot(S_E)) + P = 0.5*( L_E.dot(M).dot(L_E) + + 0.5*L_H.dot(M).dot(L_H).dot(G) + + 0.5*G.T.dot(L_H).dot(M).dot(L_H)) if self.mesh.boundary_label['LEFT'] != 'Periodic': raise ValueError("Only implemented for periodic.") @@ -269,7 +266,7 @@ def getTotalEnergy(self, G, fields): f['E'] = np.zeros(len(fields['E'])-1) f['E'][:] = fields['E'][:-1] q = self.fieldsAsStateVector(f) - return q.T.dot(V).dot(q) + return q.T.dot(P).dot(q) def reorder_by_elements(self, A): # Assumes that the original array contains all DoF ordered as: diff --git a/test/tests_fdtd/test_fdtd_1d.py b/test/tests_fdtd/test_fdtd_1d.py index 810370d..0383b34 100644 --- a/test/tests_fdtd/test_fdtd_1d.py +++ b/test/tests_fdtd/test_fdtd_1d.py @@ -243,7 +243,7 @@ def test_energy_evolution(): This requires a special operator to compute energy. ''' sp = FD1D(mesh=Mesh1D(-1.0, 1.0, 100, boundary_label="Periodic")) - dr = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=1.0) + dr = MaxwellDriver(sp, timeIntegratorType='LF2', CFL=0.7) G = dr.buildDrivedEvolutionOperator(reduceToEssentialDoF=True) s0 = 0.25 @@ -264,8 +264,8 @@ def test_energy_evolution(): # plt.plot(energyE + energyH) # plt.plot((energyH[:-1] + energyH[1:])*0.5 + energyE[:-1]) - # plt.plot(energyE) - # plt.plot(energyH) + # # plt.plot(energyE) + # # plt.plot(energyH) # plt.plot(totalEnergy) # plt.show() assert np.isclose(totalEnergy[0],totalEnergy[-1]) From e5b0c918ab808e07747d5cbe768d1208b55e3413 Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Fri, 16 Aug 2024 12:19:20 +0200 Subject: [PATCH 24/25] Minor --- maxwell/dg/dg1d.py | 20 ++++++++++++++++++++ maxwell/fd/fd1d.py | 10 ++++++---- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/maxwell/dg/dg1d.py b/maxwell/dg/dg1d.py index 500b03a..252ac96 100644 --- a/maxwell/dg/dg1d.py +++ b/maxwell/dg/dg1d.py @@ -304,6 +304,26 @@ def getEnergy(self, field): ) return energy + + def getTotalEnergy(self, G, fields): + N = self.number_of_unknowns() + NE = self.number_of_unknowns('E') + NH = self.number_of_unknowns('H') + + L_E = np.zeros((N, N)) + L_E[:NE, :NE] = np.eye(NE) + L_H = np.zeros((N, N)) + L_H[NE:, NE:] = np.eye(NH) + + M = self.buildGlobalMassMatrix() + + P = 0.5*( M + + 0.5*L_E.dot(G.T).dot(L_H).dot(M).dot(L_H) + + 0.5*L_H.dot(M).dot(L_H).dot(G).dot(L_E) + ) + + q = self.fieldsAsStateVector(fields) + return q.T.dot(P).dot(q) def buildConnectedOperators(self, element=0, neighbors=1): diff --git a/maxwell/fd/fd1d.py b/maxwell/fd/fd1d.py index abb0e43..3932734 100644 --- a/maxwell/fd/fd1d.py +++ b/maxwell/fd/fd1d.py @@ -259,12 +259,14 @@ def getTotalEnergy(self, G, fields): + 0.5*L_H.dot(M).dot(L_H).dot(G) + 0.5*G.T.dot(L_H).dot(M).dot(L_H)) - if self.mesh.boundary_label['LEFT'] != 'Periodic': + f = copy.deepcopy(fields) + + if self.mesh.boundary_label['LEFT'] == 'Periodic': + f['E'] = np.zeros(len(fields['E'])-1) + f['E'][:] = fields['E'][:-1] + else: raise ValueError("Only implemented for periodic.") - f = copy.deepcopy(fields) - f['E'] = np.zeros(len(fields['E'])-1) - f['E'][:] = fields['E'][:-1] q = self.fieldsAsStateVector(f) return q.T.dot(P).dot(q) From 63c9aa46d2bc150ddcfb843a2a883543a16eaeff Mon Sep 17 00:00:00 2001 From: Luis Manuel Diaz Angulo Date: Mon, 25 Nov 2024 08:53:41 +0100 Subject: [PATCH 25/25] Minor --- test/test_dgtd_2d.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/test_dgtd_2d.py b/test/test_dgtd_2d.py index 80e0ad6..c0cd2a3 100644 --- a/test/test_dgtd_2d.py +++ b/test/test_dgtd_2d.py @@ -18,7 +18,7 @@ def resonant_cavity_ez_field(x, y, t): return np.sin(m*np.pi*x)*np.sin(n*np.pi*y)*np.cos(w*t) def test_pec(): - N = 2 + N = 5 msh = readFromGambitFile(TEST_DATA_FOLDER + 'Maxwell2D_K146.neu') sp = Maxwell2D(N, msh, 'Centered') @@ -31,11 +31,10 @@ def test_pec(): # plt.show() for _ in range(40): - # sp.plot_field(N, driver['Ez']) - # plt.pause(0.001) - # plt.cla() - + sp.plot_field(N, driver['Ez']) + plt.pause(0.1) driver.step() + plt.cla() ez_expected = resonant_cavity_ez_field(sp.x, sp.y, driver.timeIntegrator.time) R = np.corrcoef(ez_expected, driver['Ez'])