From 4e64fc8b8775a83038ea5e8bf883b8c079c6c027 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Tue, 21 Feb 2017 12:54:00 +0000 Subject: [PATCH] Sources now use pre-calculated waveform values, instead of calculating on-the-fly. Sources now use iteration counter - no absolute time value required. Transmission line spatial step updated from using magic time step (which showed instabilities) to sqrt(3) x c x dt --- gprMax/sources.py | 139 +++++++++++++++++++++++----------------------- 1 file changed, 70 insertions(+), 69 deletions(-) diff --git a/gprMax/sources.py b/gprMax/sources.py index 62008302..3ddb6e09 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -17,6 +17,7 @@ # along with gprMax. If not, see . from copy import deepcopy +import decimal as d import numpy as np @@ -41,51 +42,67 @@ class Source(object): self.stop = None self.waveformID = None + def calculate_waveform_values(self, G): + """Calculates all waveform values for source for duration of simulation. + + Args: + G (class): Grid class instance - holds essential parameters describing the model. + """ + + self.waveformvaluesJ = np.zeros((G.iterations + 1), dtype=floattype) + self.waveformvaluesM = np.zeros((G.iterations + 1), dtype=floattype) + waveform = next(x for x in G.waveforms if x.ID == self.waveformID) + + for iteration in range(G.iterations + 1): + time = G.dt * iteration + if time >= self.start and time <= self.stop: + # Set the time of the waveform evaluation to account for any delay in the start + time -= self.start + self.waveformvaluesJ[iteration] = waveform.calculate_value(time + 0.5 * G.dt, G.dt) + self.waveformvaluesM[iteration] = waveform.calculate_value(time, G.dt) + class VoltageSource(Source): """A voltage source can be a hard source if it's resistance is zero, i.e. the time variation of the specified electric field component is prescribed. If it's resistance is non-zero it behaves as a resistive voltage source.""" def __init__(self): - super(Source, self).__init__() + super().__init__() self.resistance = None - def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): + def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G): """Updates electric field values for a voltage source. Args: - abstime (float): Absolute time. + iteration (int): Current iteration (timestep). updatecoeffsE (memory view): numpy array of electric field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. G (class): Grid class instance - holds essential parameters describing the model. """ - if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account for any delay in the start - time = abstime - self.start + if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: i = self.xcoord j = self.ycoord k = self.zcoord - waveform = next(x for x in G.waveforms if x.ID == self.waveformID) componentID = 'E' + self.polarisation if self.polarisation == 'x': if self.resistance != 0: - Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dy * G.dz)) + Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dy * G.dz)) else: - Ex[i, j, k] = -1 * waveform.amp * waveform.calculate_value(time, G.dt) / G.dx + Ex[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dx elif self.polarisation == 'y': if self.resistance != 0: - Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dx * G.dz)) + Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dx * G.dz)) else: - Ey[i, j, k] = -1 * waveform.amp * waveform.calculate_value(time, G.dt) / G.dy + Ey[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dy elif self.polarisation == 'z': if self.resistance != 0: - Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dx * G.dy)) + Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dx * G.dy)) else: - Ez[i, j, k] = -1 * waveform.amp * waveform.calculate_value(time, G.dt) / G.dz + Ez[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dz def create_material(self, G): """Create a new material at the voltage source location that adds the voltage source conductivity to the underlying parameters. @@ -106,7 +123,7 @@ class VoltageSource(Source): newmaterial.ID = material.ID + '+VoltageSource_' + str(self.resistance) newmaterial.numID = len(G.materials) newmaterial.averagable = False - newmaterial.type += ', voltage-source' + newmaterial.type += ',\nvoltage-source' # Add conductivity of voltage source to underlying conductivity if self.polarisation == 'x': @@ -124,74 +141,67 @@ class HertzianDipole(Source): """A Hertzian dipole is an additive source (electric current density).""" def __init__(self): - super(Source, self).__init__() + super().__init__() self.dl = None - def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): + def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G): """Updates electric field values for a Hertzian dipole. Args: - abstime (float): Absolute time. + iteration (int): Current iteration (timestep). updatecoeffsE (memory view): numpy array of electric field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. G (class): Grid class instance - holds essential parameters describing the model. """ - if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account for any delay in the start - time = abstime - self.start + if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: i = self.xcoord j = self.ycoord k = self.zcoord - waveform = next(x for x in G.waveforms if x.ID == self.waveformID) componentID = 'E' + self.polarisation if self.polarisation == 'x': - Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) elif self.polarisation == 'y': - Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) elif self.polarisation == 'z': - Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) class MagneticDipole(Source): """A magnetic dipole is an additive source (magnetic current density).""" def __init__(self): - super(Source, self).__init__() + super().__init__() - def update_magnetic(self, abstime, updatecoeffsH, ID, Hx, Hy, Hz, G): + def update_magnetic(self, iteration, updatecoeffsH, ID, Hx, Hy, Hz, G): """Updates magnetic field values for a magnetic dipole. Args: - abstime (float): Absolute time. + iteration (int): Current iteration (timestep). updatecoeffsH (memory view): numpy array of magnetic field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Hx, Hy, Hz (memory view): numpy array of magnetic field values. G (class): Grid class instance - holds essential parameters describing the model. """ - if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account for any delay in the start - time = abstime - self.start + if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: i = self.xcoord j = self.ycoord k = self.zcoord - waveform = next(x for x in G.waveforms if x.ID == self.waveformID) componentID = 'H' + self.polarisation if self.polarisation == 'x': - Hx[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy * G.dz)) + Hx[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) elif self.polarisation == 'y': - - Hy[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy * G.dz)) + Hy[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) elif self.polarisation == 'z': - Hz[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy * G.dz)) + Hz[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) class TransmissionLine(Source): @@ -203,15 +213,15 @@ class TransmissionLine(Source): G (class): Grid class instance - holds essential parameters describing the model. """ - super(Source, self).__init__() + super().__init__() self.resistance = None # Coefficients for ABC termination of end of the transmission line self.abcv0 = 0 self.abcv1 = 0 - # Spatial step of transmission line (based on magic time step for dispersionless behaviour) - self.dl = c * G.dt + # Spatial step of transmission line (N.B if the magic time step is used it results in instabilities for certain impedances) + self.dl = np.sqrt(3) * c * G.dt # Number of cells in the transmission line (initially a long line to calculate incident voltage and current); consider putting ABCs/PML at end self.nl = round_value(0.667 * G.iterations) @@ -236,14 +246,11 @@ class TransmissionLine(Source): G (class): Grid class instance - holds essential parameters describing the model. """ - abstime = 0 - for timestep in range(G.iterations): - self.Vinc[timestep] = self.voltage[self.antpos] - self.Iinc[timestep] = self.current[self.antpos] - self.update_voltage(abstime, G) - abstime += 0.5 * G.dt - self.update_current(abstime, G) - abstime += 0.5 * G.dt + for iteration in range(G.iterations): + self.Iinc[iteration] = self.current[self.antpos] + self.Vinc[iteration] = self.voltage[self.antpos] + self.update_current(iteration, G) + self.update_voltage(iteration, G) # Shorten number of cells in the transmission line before use with main grid self.nl = self.antpos + 1 @@ -261,11 +268,11 @@ class TransmissionLine(Source): self.abcv0 = self.voltage[0] self.abcv1 = self.voltage[1] - def update_voltage(self, time, G): + def update_voltage(self, iteration, G): """Updates voltage values along the transmission line. Args: - time (float): Absolute time. + iteration (int): Current iteration (timestep). G (class): Grid class instance - holds essential parameters describing the model. """ @@ -273,17 +280,16 @@ class TransmissionLine(Source): self.voltage[1:self.nl] -= self.resistance * (c * G.dt / self.dl) * (self.current[1:self.nl] - self.current[0:self.nl - 1]) # Update the voltage at the position of the one-way injector excitation - waveform = next(x for x in G.waveforms if x.ID == self.waveformID) - self.voltage[self.srcpos] += (c * G.dt / self.dl) * waveform.amp * waveform.calculate_value(time - 0.5 * G.dt, G.dt) + self.voltage[self.srcpos] += (c * G.dt / self.dl) * self.waveformvaluesJ[iteration] # Update ABC before updating current self.update_abc(G) - def update_current(self, time, G): + def update_current(self, iteration, G): """Updates current values along the transmission line. Args: - time (float): Absolute time. + iteration (int): Current iteration (timestep). G (class): Grid class instance - holds essential parameters describing the model. """ @@ -291,28 +297,25 @@ class TransmissionLine(Source): self.current[0:self.nl - 1] -= (1 / self.resistance) * (c * G.dt / self.dl) * (self.voltage[1:self.nl] - self.voltage[0:self.nl - 1]) # Update the current one cell before the position of the one-way injector excitation - waveform = next(x for x in G.waveforms if x.ID == self.waveformID) - self.current[self.srcpos - 1] += (c * G.dt / self.dl) * waveform.amp * waveform.calculate_value(time - 0.5 * G.dt, G.dt) * (1 / self.resistance) + self.current[self.srcpos - 1] += (1 / self.resistance) * (c * G.dt / self.dl) * self.waveformvaluesM[iteration] - def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): + def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G): """Updates electric field value in the main grid from voltage value in the transmission line. Args: - abstime (float): Absolute time. + iteration (int): Current iteration (timestep). updatecoeffsE (memory view): numpy array of electric field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. G (class): Grid class instance - holds essential parameters describing the model. """ - if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account for any delay in the start - time = abstime - self.start + if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: i = self.xcoord j = self.ycoord k = self.zcoord - self.update_voltage(time, G) + self.update_voltage(iteration, G) if self.polarisation == 'x': Ex[i, j, k] = - self.voltage[self.antpos] / G.dx @@ -323,34 +326,32 @@ class TransmissionLine(Source): elif self.polarisation == 'z': Ez[i, j, k] = - self.voltage[self.antpos] / G.dz - def update_magnetic(self, abstime, updatecoeffsH, ID, Hx, Hy, Hz, G): + def update_magnetic(self, iteration, updatecoeffsH, ID, Hx, Hy, Hz, G): """Updates current value in transmission line from magnetic field values in the main grid. Args: - abstime (float): Absolute time. + iteration (int): Current iteration (timestep). updatecoeffsH (memory view): numpy array of magnetic field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Hx, Hy, Hz (memory view): numpy array of magnetic field values. G (class): Grid class instance - holds essential parameters describing the model. """ - if abstime >= self.start and abstime <= self.stop: - # Set the time of the waveform evaluation to account for any delay in the start - time = abstime - self.start + if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: i = self.xcoord j = self.ycoord k = self.zcoord if self.polarisation == 'x': - self.current[self.antpos] = Ix(i, j, k, G.Hy, G.Hz, G) + self.current[self.antpos] = Ix(i, j, k, G.Hx, G.Hy, G.Hz, G) elif self.polarisation == 'y': - self.current[self.antpos] = Iy(i, j, k, G.Hx, G.Hz, G) + self.current[self.antpos] = Iy(i, j, k, G.Hx, G.Hy, G.Hz, G) elif self.polarisation == 'z': - self.current[self.antpos] = Iz(i, j, k, G.Hx, G.Hy, G) + self.current[self.antpos] = Iz(i, j, k, G.Hx, G.Hy, G.Hz, G) - self.update_current(time, G) + self.update_current(iteration, G) class PlaneWave(Source):