multi cmds to functions

这个提交包含在:
jasminium
2016-06-09 18:05:42 +01:00
父节点 39fe2675a3
当前提交 92bf59bc36
共有 5 个文件被更改,包括 928 次插入824 次删除

查看文件

@@ -68,7 +68,7 @@ def main():
def api(inputfile, n=1, mpi=False, benchmark=False, geometry_only=False, geometry_fixed=False, write_processed=False, opt_taguchi=False): def api(inputfile, n=1, mpi=False, benchmark=False, geometry_only=False, geometry_fixed=False, write_processed=False, opt_taguchi=False):
"""If installed as a module this is the entry point.""" """If installed as a module this is the entry point."""
class ImportArguments: class ImportArguments:
pass pass
@@ -88,11 +88,11 @@ def api(inputfile, n=1, mpi=False, benchmark=False, geometry_only=False, geometr
def run_main(args): def run_main(args):
"""Top-level function that controls what mode of simulation (standard/optimsation/benchmark etc...) is run. """Top-level function that controls what mode of simulation (standard/optimsation/benchmark etc...) is run.
Args: Args:
args (dict): Namespace with input arguments from command line or api. args (dict): Namespace with input arguments from command line or api.
""" """
numbermodelruns = args.n numbermodelruns = args.n
inputdirectory = os.path.dirname(os.path.abspath(args.inputfile)) inputdirectory = os.path.dirname(os.path.abspath(args.inputfile))
inputfile = os.path.abspath(os.path.join(inputdirectory, os.path.basename(args.inputfile))) inputfile = os.path.abspath(os.path.join(inputdirectory, os.path.basename(args.inputfile)))

查看文件

@@ -25,32 +25,32 @@ from gprMax.utilities import ListStream
def process_python_include_code(inputfile, usernamespace): def process_python_include_code(inputfile, usernamespace):
"""Looks for and processes any Python code found in the input file. It will ignore any lines that are comments, i.e. begin with a double hash (##), and any blank lines. It will also ignore any lines that do not begin with a hash (#) after it has processed Python commands. It will also process any include commands and insert the contents of the included file at that location. """Looks for and processes any Python code found in the input file. It will ignore any lines that are comments, i.e. begin with a double hash (##), and any blank lines. It will also ignore any lines that do not begin with a hash (#) after it has processed Python commands. It will also process any include commands and insert the contents of the included file at that location.
Args: Args:
inputfile (str): Name of the input file to open. inputfile (str): Name of the input file to open.
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
Returns: Returns:
processedlines (list): Input commands after Python processing. processedlines (list): Input commands after Python processing.
""" """
userpython = False userpython = False
with open(inputfile, 'r') as f: with open(inputfile, 'r') as f:
# Strip out any newline characters and comments that must begin with double hashes # Strip out any newline characters and comments that must begin with double hashes
inputlines = [line.rstrip() for line in f if(not line.startswith('##') and line.rstrip('\n'))] inputlines = [line.rstrip() for line in f if(not line.startswith('##') and line.rstrip('\n'))]
# List to hold final processed commands # List to hold final processed commands
processedlines = [] processedlines = []
x = 0 x = 0
while(x < len(inputlines)): while(x < len(inputlines)):
# Process any Python code # Process any Python code
if(inputlines[x].startswith('#python:')): if(inputlines[x].startswith('#python:')):
# Save stdout location to restore later # Save stdout location to restore later
stdout = sys.stdout stdout = sys.stdout
# String to hold Python code to be executed # String to hold Python code to be executed
pythoncode = '' pythoncode = ''
x += 1 x += 1
@@ -66,36 +66,36 @@ def process_python_include_code(inputfile, usernamespace):
sys.stdout = codeout = ListStream() sys.stdout = codeout = ListStream()
# Execute code block & make available only usernamespace # Execute code block & make available only usernamespace
exec(pythoncompiledcode, usernamespace) exec(pythoncompiledcode, usernamespace)
# Now strip out any lines that don't begin with a hash command # Now strip out any lines that don't begin with a hash command
codeproc = [line + ('\n') for line in codeout.data if(line.startswith('#'))] codeproc = [line + ('\n') for line in codeout.data if(line.startswith('#'))]
# Add processed Python code to list # Add processed Python code to list
processedlines.extend(codeproc) processedlines.extend(codeproc)
# Reset stdio # Reset stdio
sys.stdout = stdout sys.stdout = stdout
# Process any include commands # Process any include commands
elif(inputlines[x].startswith('#include_file:')): elif(inputlines[x].startswith('#include_file:')):
includefile = inputlines[x].split() includefile = inputlines[x].split()
if len(includefile) != 2: if len(includefile) != 2:
raise CmdInputError('#include_file requires exactly one parameter') raise CmdInputError('#include_file requires exactly one parameter')
includefile = includefile[1] includefile = includefile[1]
# See if file exists at specified path and if not try input file directory # See if file exists at specified path and if not try input file directory
if not os.path.isfile(includefile): if not os.path.isfile(includefile):
includefile = os.path.join(usernamespace['inputdirectory'], includefile) includefile = os.path.join(usernamespace['inputdirectory'], includefile)
with open(includefile, 'r') as f: with open(includefile, 'r') as f:
# Strip out any newline characters and comments that must begin with double hashes # Strip out any newline characters and comments that must begin with double hashes
includelines = [includeline.rstrip() + '\n' for includeline in f if(not includeline.startswith('##') and includeline.rstrip('\n'))] includelines = [includeline.rstrip() + '\n' for includeline in f if(not includeline.startswith('##') and includeline.rstrip('\n'))]
# Add lines from include file to list # Add lines from include file to list
processedlines.extend(includelines) processedlines.extend(includelines)
# Add any other commands to list # Add any other commands to list
elif(inputlines[x].startswith('#')): elif(inputlines[x].startswith('#')):
# Add gprMax command to list # Add gprMax command to list
@@ -103,20 +103,20 @@ def process_python_include_code(inputfile, usernamespace):
processedlines.append(inputlines[x]) processedlines.append(inputlines[x])
x += 1 x += 1
return processedlines return processedlines
def write_processed_file(inputfile, modelrun, numbermodelruns, processedlines): def write_processed_file(inputfile, modelrun, numbermodelruns, processedlines):
"""Writes an input file after any Python code and include commands in the original input file have been processed. """Writes an input file after any Python code and include commands in the original input file have been processed.
Args: Args:
inputfile (str): Name of the input file to open. inputfile (str): Name of the input file to open.
modelrun (int): Current model run number. modelrun (int): Current model run number.
numbermodelruns (int): Total number of model runs. numbermodelruns (int): Total number of model runs.
processedlines (list): Input commands after after processing any Python code and include commands. processedlines (list): Input commands after after processing any Python code and include commands.
""" """
if numbermodelruns == 1: if numbermodelruns == 1:
processedfile = os.path.splitext(inputfile)[0] + '_processed.in' processedfile = os.path.splitext(inputfile)[0] + '_processed.in'
else: else:
@@ -131,11 +131,11 @@ def write_processed_file(inputfile, modelrun, numbermodelruns, processedlines):
def check_cmd_names(processedlines, checkessential=True): def check_cmd_names(processedlines, checkessential=True):
"""Checks the validity of commands, i.e. are they gprMax commands, and that all essential commands are present. """Checks the validity of commands, i.e. are they gprMax commands, and that all essential commands are present.
Args: Args:
processedlines (list): Input commands after Python processing. processedlines (list): Input commands after Python processing.
checkessential (boolean): Perform check to see that all essential commands are present. checkessential (boolean): Perform check to see that all essential commands are present.
Returns: Returns:
singlecmds (dict): Commands that can only occur once in the model. singlecmds (dict): Commands that can only occur once in the model.
multiplecmds (dict): Commands that can have multiple instances in the model. multiplecmds (dict): Commands that can have multiple instances in the model.
@@ -150,13 +150,13 @@ def check_cmd_names(processedlines, checkessential=True):
singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi'], 'None') singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi'], 'None')
# Commands that there can be multiple instances of in a model - these will be lists within the dictionary # Commands that there can be multiple instances of in a model - these will be lists within the dictionary
multiplecmds = {key: [] for key in ['#geometry_view', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#rx', '#rx_box', '#snapshot', '#pml_cfs']} multiplecmds = {key: [] for key in ['#geometry_view', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#tem_transmission_line', '#rx', '#rx_box', '#snapshot', '#pml_cfs']}
# Geometry object building commands that there can be multiple instances of in a model - these will be lists within the dictionary # Geometry object building commands that there can be multiple instances of in a model - these will be lists within the dictionary
geometrycmds = ['#geometry_objects_file', '#edge', '#plate', '#triangle', '#box', '#sphere', '#cylinder', '#cylindrical_sector', '#fractal_box', '#add_surface_roughness', '#add_surface_water', '#add_grass'] geometrycmds = ['#geometry_objects_file', '#edge', '#plate', '#triangle', '#box', '#sphere', '#cylinder', '#cylindrical_sector', '#fractal_box', '#add_surface_roughness', '#add_surface_water', '#add_grass']
# List to store all geometry object commands in order from input file # List to store all geometry object commands in order from input file
geometry = [] geometry = []
# Check if command names are valid, if essential commands are present, and add command parameters to appropriate dictionary values or lists # Check if command names are valid, if essential commands are present, and add command parameters to appropriate dictionary values or lists
countessentialcmds = 0 countessentialcmds = 0
lindex = 0 lindex = 0
@@ -171,22 +171,22 @@ def check_cmd_names(processedlines, checkessential=True):
# Count essential commands # Count essential commands
if cmdname in essentialcmds: if cmdname in essentialcmds:
countessentialcmds += 1 countessentialcmds += 1
# Assign command parameters as values to dictionary keys # Assign command parameters as values to dictionary keys
if cmdname in singlecmds: if cmdname in singlecmds:
if singlecmds[cmdname] == 'None': if singlecmds[cmdname] == 'None':
singlecmds[cmdname] = cmd[1].strip(' \t\n') singlecmds[cmdname] = cmd[1].strip(' \t\n')
else: else:
raise CmdInputError('You can only have one ' + cmdname + ' commmand in your model') raise CmdInputError('You can only have one ' + cmdname + ' commmand in your model')
elif cmdname in multiplecmds: elif cmdname in multiplecmds:
multiplecmds[cmdname].append(cmd[1].strip(' \t\n')) multiplecmds[cmdname].append(cmd[1].strip(' \t\n'))
elif cmdname in geometrycmds: elif cmdname in geometrycmds:
geometry.append(processedlines[lindex].strip(' \t\n')) geometry.append(processedlines[lindex].strip(' \t\n'))
lindex += 1 lindex += 1
if checkessential: if checkessential:
if (countessentialcmds < len(essentialcmds)): if (countessentialcmds < len(essentialcmds)):
raise CmdInputError('Your input file is missing essential gprMax commands required to run a model. Essential commands are: ' + ', '.join(essentialcmds)) raise CmdInputError('Your input file is missing essential gprMax commands required to run a model. Essential commands are: ' + ', '.join(essentialcmds))

文件差异内容过多而无法显示 加载差异

查看文件

@@ -20,28 +20,34 @@ import numpy as np
from gprMax.constants import e0, m0, floattype, complextype from gprMax.constants import e0, m0, floattype, complextype
# Look up for material ids
material_ids = {
'pec': 0,
'freespace': 1
}
class Material(object): class Material(object):
"""Materials, their properties and update coefficients.""" """Materials, their properties and update coefficients."""
# Maximum number of dispersive material poles in a model # Maximum number of dispersive material poles in a model
maxpoles = 0 maxpoles = 0
# Types of material # Types of material
types = ['standard', 'debye', 'lorentz', 'drude'] types = ['standard', 'debye', 'lorentz', 'drude']
# Properties of water from: http://dx.doi.org/10.1109/TGRS.2006.873208 # Properties of water from: http://dx.doi.org/10.1109/TGRS.2006.873208
waterer = 80.1 waterer = 80.1
watereri = 4.9 watereri = 4.9
waterdeltaer = waterer - watereri waterdeltaer = waterer - watereri
watertau = 9.231e-12 watertau = 9.231e-12
# Properties of grass from: http://dx.doi.org/10.1007/BF00902994 # Properties of grass from: http://dx.doi.org/10.1007/BF00902994
grasser = 18.5087 grasser = 18.5087
grasseri = 12.7174 grasseri = 12.7174
grassdeltaer = grasser - grasseri grassdeltaer = grasser - grasseri
grasstau = 1.0793e-11 grasstau = 1.0793e-11
def __init__(self, numID, ID, G): def __init__(self, numID, ID, G):
""" """
Args: Args:
@@ -49,7 +55,7 @@ class Material(object):
ID (str): Name of the material. ID (str): Name of the material.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
self.numID = numID self.numID = numID
self.ID = ID self.ID = ID
self.type = 'standard' self.type = 'standard'
@@ -61,20 +67,20 @@ class Material(object):
self.se = 0.0 self.se = 0.0
self.mr = 1.0 self.mr = 1.0
self.sm = 0.0 self.sm = 0.0
# Parameters for dispersive materials # Parameters for dispersive materials
self.poles = 0 self.poles = 0
self.deltaer = [] self.deltaer = []
self.tau = [] self.tau = []
self.alpha = [] self.alpha = []
def calculate_update_coeffsH(self, G): def calculate_update_coeffsH(self, G):
"""Calculates the magnetic update coefficients of the material. """Calculates the magnetic update coefficients of the material.
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
HA = (m0*self.mr / G.dt) + 0.5*self.sm HA = (m0*self.mr / G.dt) + 0.5*self.sm
HB = (m0*self.mr / G.dt) - 0.5*self.sm HB = (m0*self.mr / G.dt) - 0.5*self.sm
self.DA = HB / HA self.DA = HB / HA
@@ -86,11 +92,11 @@ class Material(object):
# Calculate electric update coefficients # Calculate electric update coefficients
def calculate_update_coeffsE(self, G): def calculate_update_coeffsE(self, G):
"""Calculates the electric update coefficients of the material. """Calculates the electric update coefficients of the material.
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
# The implementation of the dispersive material modelling comes from the derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549 # The implementation of the dispersive material modelling comes from the derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549
if self.maxpoles > 0: if self.maxpoles > 0:
self.w = np.zeros(self.maxpoles, dtype=complextype) self.w = np.zeros(self.maxpoles, dtype=complextype)
@@ -99,7 +105,7 @@ class Material(object):
self.zt2 = np.zeros(self.maxpoles, dtype=complextype) self.zt2 = np.zeros(self.maxpoles, dtype=complextype)
self.eqt = np.zeros(self.maxpoles, dtype=complextype) self.eqt = np.zeros(self.maxpoles, dtype=complextype)
self.eqt2 = np.zeros(self.maxpoles, dtype=complextype) self.eqt2 = np.zeros(self.maxpoles, dtype=complextype)
for x in range(self.poles): for x in range(self.poles):
if self.type == 'debye': if self.type == 'debye':
self.w[x] = self.deltaer[x] / self.tau[x] self.w[x] = self.deltaer[x] / self.tau[x]
@@ -117,7 +123,7 @@ class Material(object):
self.se += wp2 / self.alpha[x] self.se += wp2 / self.alpha[x]
self.w[x] = - (wp2 / self.alpha[x]) self.w[x] = - (wp2 / self.alpha[x])
self.q[x] = - self.alpha[x] self.q[x] = - self.alpha[x]
self.eqt[x] = np.exp(self.q[x] * G.dt) self.eqt[x] = np.exp(self.q[x] * G.dt)
self.eqt2[x] = np.exp(self.q[x] * (G.dt / 2)) self.eqt2[x] = np.exp(self.q[x] * (G.dt / 2))
self.zt[x] = (self.w[x] / self.q[x]) * (1 - self.eqt[x]) / G.dt self.zt[x] = (self.w[x] / self.q[x]) * (1 - self.eqt[x]) / G.dt
@@ -125,7 +131,7 @@ class Material(object):
EA = (e0*self.er / G.dt) + 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real) EA = (e0*self.er / G.dt) + 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real)
EB = (e0*self.er / G.dt) - 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real) EB = (e0*self.er / G.dt) - 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real)
else: else:
EA = (e0*self.er / G.dt) + 0.5*self.se EA = (e0*self.er / G.dt) + 0.5*self.se
EB = (e0*self.er / G.dt) - 0.5*self.se EB = (e0*self.er / G.dt) - 0.5*self.se
@@ -146,7 +152,7 @@ class Material(object):
class PeplinskiSoil(object): class PeplinskiSoil(object):
"""Soil objects that are characterised according to a mixing model by Peplinski (http://dx.doi.org/10.1109/36.387598).""" """Soil objects that are characterised according to a mixing model by Peplinski (http://dx.doi.org/10.1109/36.387598)."""
def __init__(self, ID, sandfraction, clayfraction, bulkdensity, sandpartdensity, watervolfraction): def __init__(self, ID, sandfraction, clayfraction, bulkdensity, sandpartdensity, watervolfraction):
""" """
Args: Args:
@@ -157,7 +163,7 @@ class PeplinskiSoil(object):
sandpartdensity (float): Density of the sand particles in the soil (g/cm3). sandpartdensity (float): Density of the sand particles in the soil (g/cm3).
watervolfraction (float): Two numbers that specify a range for the volumetric water fraction of the soil. watervolfraction (float): Two numbers that specify a range for the volumetric water fraction of the soil.
""" """
self.ID = ID self.ID = ID
self.S = sandfraction self.S = sandfraction
self.C = clayfraction self.C = clayfraction
@@ -168,28 +174,28 @@ class PeplinskiSoil(object):
def calculate_debye_properties(self, nbins, G): def calculate_debye_properties(self, nbins, G):
"""Calculates the real and imaginery part of a Debye model for the soil as well as a conductivity. It uses a semi-empirical model (http://dx.doi.org/10.1109/36.387598). """Calculates the real and imaginery part of a Debye model for the soil as well as a conductivity. It uses a semi-empirical model (http://dx.doi.org/10.1109/36.387598).
Args: Args:
nbins (int): Number of bins to use to create the different materials. nbins (int): Number of bins to use to create the different materials.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
# Debye model properties of water # Debye model properties of water
f = 1.3e9 f = 1.3e9
w = 2 * np.pi * f w = 2 * np.pi * f
erealw = Material.watereri + ((Material.waterdeltaer) / (1 + (w * Material.watertau)**2)) erealw = Material.watereri + ((Material.waterdeltaer) / (1 + (w * Material.watertau)**2))
eimagw = w * Material.watertau * ((Material.waterdeltaer) / (1 + (w * Material.watertau)**2)) eimagw = w * Material.watertau * ((Material.waterdeltaer) / (1 + (w * Material.watertau)**2))
a = 0.65 # Experimentally derived constant a = 0.65 # Experimentally derived constant
es = (1.01 + 0.44 * self.rs)**2 - 0.062 es = (1.01 + 0.44 * self.rs)**2 - 0.062
b1 = 1.2748 - 0.519 * self.S - 0.152 * self.C b1 = 1.2748 - 0.519 * self.S - 0.152 * self.C
b2 = 1.33797 - 0.603 * self.S - 0.166 * self.C b2 = 1.33797 - 0.603 * self.S - 0.166 * self.C
# For frequencies in the range 0.3GHz to 1.3GHz # For frequencies in the range 0.3GHz to 1.3GHz
sigf1 = 0.0467 + 0.2204 * self.rb - 0.411 * self.S + 0.6614 * self.C sigf1 = 0.0467 + 0.2204 * self.rb - 0.411 * self.S + 0.6614 * self.C
# For frequencies in the range 1.4GHz to 18GHz # For frequencies in the range 1.4GHz to 18GHz
sigf2 = -1.645 + 1.939 * self.rb - 2.25622 * self.S + 1.594 * self.C sigf2 = -1.645 + 1.939 * self.rb - 2.25622 * self.S + 1.594 * self.C
# Generate a set of bins based on the given volumetric water fraction values # Generate a set of bins based on the given volumetric water fraction values
mubins = np.linspace(self.mu[0], self.mu[1], nbins + 1) mubins = np.linspace(self.mu[0], self.mu[1], nbins + 1)
# Generate a range of volumetric water fraction values the mid-point of each bin to make materials from # Generate a range of volumetric water fraction values the mid-point of each bin to make materials from
@@ -202,13 +208,13 @@ class PeplinskiSoil(object):
er1 = (1 + (self.rb/self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1/a) er1 = (1 + (self.rb/self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1/a)
# Real part for frequencies in the range 0.3GHz to 1.3GHz # Real part for frequencies in the range 0.3GHz to 1.3GHz
er2 = 1.15 * er1 - 0.68 er2 = 1.15 * er1 - 0.68
# Imaginary part for frequencies in the range 0.3GHz to 1.3GHz # Imaginary part for frequencies in the range 0.3GHz to 1.3GHz
eri = er2 - (muiter[0]**(b2/a) * Material.waterdeltaer) eri = er2 - (muiter[0]**(b2/a) * Material.waterdeltaer)
# Effective conductivity # Effective conductivity
sig = muiter[0]**(b2/a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0])) sig = muiter[0]**(b2/a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0]))
# Check to see if the material already exists before creating a new one # Check to see if the material already exists before creating a new one
requiredID = '|{:.4f}|'.format(float(muiter[0])) requiredID = '|{:.4f}|'.format(float(muiter[0]))
material = next((x for x in G.materials if x.ID == requiredID), None) material = next((x for x in G.materials if x.ID == requiredID), None)
@@ -225,7 +231,7 @@ class PeplinskiSoil(object):
m.deltaer.append(er2 - m.er) m.deltaer.append(er2 - m.er)
m.tau.append(Material.watertau) m.tau.append(Material.watertau)
G.materials.append(m) G.materials.append(m)
muiter.iternext() muiter.iternext()

查看文件

@@ -20,9 +20,11 @@ from copy import deepcopy
import numpy as np import numpy as np
from gprMax.constants import c, floattype from .constants import c, floattype
from gprMax.grid import Ix, Iy, Iz from .grid import Ix, Iy, Iz
from gprMax.utilities import round_value from .utilities import round_value
from .materials import material_ids
from .exceptions import GeneralError
class Source(object): class Source(object):
@@ -43,15 +45,18 @@ class Source(object):
class VoltageSource(Source): class VoltageSource(Source):
"""The 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.""" """The 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): def __init__(self):
super(Source, self).__init__() super().__init__()
self.resistance = None self.resistance = None
def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G):
"""Updates electric field values for a voltage source. """Updates electric field values for a voltage source.
Args: Args:
abstime (float): Absolute time. abstime (float): Absolute time.
updatecoeffsE (memory view): numpy array of electric field update coefficients. updatecoeffsE (memory view): numpy array of electric field update coefficients.
@@ -59,7 +64,7 @@ class VoltageSource(Source):
Ex, Ey, Ez (memory view): numpy array of electric field values. Ex, Ey, Ez (memory view): numpy array of electric field values.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
if abstime >= self.start and abstime <= self.stop: if abstime >= self.start and abstime <= self.stop:
# Set the time of the waveform evaluation to account for any delay in the start # Set the time of the waveform evaluation to account for any delay in the start
time = abstime - self.start time = abstime - self.start
@@ -67,7 +72,7 @@ class VoltageSource(Source):
j = self.ycoord j = self.ycoord
k = self.zcoord k = self.zcoord
waveform = next(x for x in G.waveforms if x.ID == self.waveformID) waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
if self.polarisation is 'x': if self.polarisation is 'x':
if self.resistance != 0: if self.resistance != 0:
componentID = 'E' + self.polarisation componentID = 'E' + self.polarisation
@@ -91,7 +96,7 @@ class VoltageSource(Source):
def create_material(self, G): def create_material(self, G):
"""Create a new material at the voltage source location that adds the voltage source conductivity to the underlying parameters. """Create a new material at the voltage source location that adds the voltage source conductivity to the underlying parameters.
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
@@ -100,7 +105,7 @@ class VoltageSource(Source):
i = self.xcoord i = self.xcoord
j = self.ycoord j = self.ycoord
k = self.zcoord k = self.zcoord
componentID = 'E' + self.polarisation componentID = 'E' + self.polarisation
requirednumID = G.ID[G.IDlookup[componentID], i, j, k] requirednumID = G.ID[G.IDlookup[componentID], i, j, k]
material = next(x for x in G.materials if x.numID == requirednumID) material = next(x for x in G.materials if x.numID == requirednumID)
@@ -108,7 +113,7 @@ class VoltageSource(Source):
newmaterial.ID = material.ID + '+VoltageSource_' + str(self.resistance) newmaterial.ID = material.ID + '+VoltageSource_' + str(self.resistance)
newmaterial.numID = len(G.materials) newmaterial.numID = len(G.materials)
newmaterial.average = False newmaterial.average = False
# Add conductivity of voltage source to underlying conductivity # Add conductivity of voltage source to underlying conductivity
if self.polarisation == 'x': if self.polarisation == 'x':
newmaterial.se += G.dx / (self.resistance * G.dy * G.dz) newmaterial.se += G.dx / (self.resistance * G.dy * G.dz)
@@ -123,13 +128,13 @@ class VoltageSource(Source):
class HertzianDipole(Source): class HertzianDipole(Source):
"""The Hertzian dipole is an additive source (electric current density).""" """The Hertzian dipole is an additive source (electric current density)."""
def __init__(self): def __init__(self):
super(Source, self).__init__() super().__init__()
def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G): def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G):
"""Updates electric field values for a Hertzian dipole. """Updates electric field values for a Hertzian dipole.
Args: Args:
abstime (float): Absolute time. abstime (float): Absolute time.
updatecoeffsE (memory view): numpy array of electric field update coefficients. updatecoeffsE (memory view): numpy array of electric field update coefficients.
@@ -137,7 +142,7 @@ class HertzianDipole(Source):
Ex, Ey, Ez (memory view): numpy array of electric field values. Ex, Ey, Ez (memory view): numpy array of electric field values.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
if abstime >= self.start and abstime <= self.stop: if abstime >= self.start and abstime <= self.stop:
# Set the time of the waveform evaluation to account for any delay in the start # Set the time of the waveform evaluation to account for any delay in the start
time = abstime - self.start time = abstime - self.start
@@ -145,7 +150,7 @@ class HertzianDipole(Source):
j = self.ycoord j = self.ycoord
k = self.zcoord k = self.zcoord
waveform = next(x for x in G.waveforms if x.ID == self.waveformID) waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
if self.polarisation is 'x': if self.polarisation is 'x':
componentID = 'E' + self.polarisation componentID = 'E' + self.polarisation
Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dy * G.dz)) Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dy * G.dz))
@@ -161,13 +166,13 @@ class HertzianDipole(Source):
class MagneticDipole(Source): class MagneticDipole(Source):
"""The magnetic dipole is an additive source (magnetic current density).""" """The magnetic dipole is an additive source (magnetic current density)."""
def __init__(self): def __init__(self):
super(Source, self).__init__() super().__init__()
def update_magnetic(self, abstime, updatecoeffsH, ID, Hx, Hy, Hz, G): def update_magnetic(self, abstime, updatecoeffsH, ID, Hx, Hy, Hz, G):
"""Updates magnetic field values for a magnetic dipole. """Updates magnetic field values for a magnetic dipole.
Args: Args:
abstime (float): Absolute time. abstime (float): Absolute time.
updatecoeffsH (memory view): numpy array of magnetic field update coefficients. updatecoeffsH (memory view): numpy array of magnetic field update coefficients.
@@ -175,7 +180,7 @@ class MagneticDipole(Source):
Hx, Hy, Hz (memory view): numpy array of magnetic field values. Hx, Hy, Hz (memory view): numpy array of magnetic field values.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
if abstime >= self.start and abstime <= self.stop: if abstime >= self.start and abstime <= self.stop:
# Set the time of the waveform evaluation to account for any delay in the start # Set the time of the waveform evaluation to account for any delay in the start
time = abstime - self.start time = abstime - self.start
@@ -183,7 +188,7 @@ class MagneticDipole(Source):
j = self.ycoord j = self.ycoord
k = self.zcoord k = self.zcoord
waveform = next(x for x in G.waveforms if x.ID == self.waveformID) waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
if self.polarisation is 'x': if self.polarisation is 'x':
Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
@@ -195,45 +200,60 @@ class MagneticDipole(Source):
class TransmissionLine(Source): class TransmissionLine(Source):
"""The transmission line source is a one-dimensional transmission line which is attached virtually to a grid cell.""" """The transmission line source is a one-dimensional transmission
line which is attached virtually to a grid cell.
"""
def __init__(self, G): def __init__(self, G):
""" """
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
super(Source, self).__init__() super().__init__()
self.resistance = None self.resistance = None
# Coefficients for ABC termination of end of the transmission line # Coefficients for ABC termination of end of the transmission line
self.abcv0 = 0 self.abcv0 = 0
self.abcv1 = 0 self.abcv1 = 0
# Spatial step of transmission line (based on magic time step for dispersionless behaviour) # Spatial step of transmission line (based on magic time step for
# dispersionless behaviour)
self.dl = c * G.dt self.dl = c * G.dt
# Number of nodes in the transmission line (initially a long line to calculate incident voltage and current); consider putting ABCs/PML at end # Number of nodes 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) self.nl = round_value(0.667 * G.iterations)
# Nodal position of the one-way injector excitation in the transmission line # Nodal position of the one-way injector excitation in the
# transmission line
self.srcpos = 5 self.srcpos = 5
# Nodal position of where line connects to antenna/main grid # Nodal position of where line connects to antenna/main grid
self.antpos = 10 self.antpos = 10
self.voltage = np.zeros(self.nl, dtype=floattype) self.voltage = np.zeros(self.nl, dtype=floattype)
self.current = np.zeros(self.nl, dtype=floattype) self.current = np.zeros(self.nl, dtype=floattype)
self.Vinc = np.zeros(G.iterations, dtype=floattype) self.Vinc = np.zeros(G.iterations, dtype=floattype)
self.Iinc = np.zeros(G.iterations, dtype=floattype) self.Iinc = np.zeros(G.iterations, dtype=floattype)
self.type_str = 'TransmissionLine'
def setID(self):
self.ID = self.type_str + '(' + str(self.xcoord) + ',' + str(self.ycoord) + ',' + str(self.zcoord) + ')'
def calculate_incident_V_I(self, G): def calculate_incident_V_I(self, G):
"""Calculates the incident voltage and current with a long length transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415 """Calculates the incident voltage and current with a long length
transmission line not connected to the main grid from:
http://dx.doi.org/10.1002/mop.10415
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
abstime = 0 abstime = 0
for timestep in range(G.iterations): for timestep in range(G.iterations):
self.Vinc[timestep] = self.voltage[self.antpos] self.Vinc[timestep] = self.voltage[self.antpos]
@@ -243,33 +263,36 @@ class TransmissionLine(Source):
self.update_current(abstime, G) self.update_current(abstime, G)
abstime += 0.5 * G.dt abstime += 0.5 * G.dt
# Shorten number of nodes in the transmission line before use with main grid # Shorten number of nodes in the transmission line before use
# with main grid
self.nl = self.antpos + 1 self.nl = self.antpos + 1
def update_abc(self, G): def update_abc(self, G):
"""Updates absorbing boundary condition at end of the transmission line. """Updates absorbing boundary condition at end of the transmission line.
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
h = (c * G.dt - self.dl) / (c * G.dt + self.dl) h = (c * G.dt - self.dl) / (c * G.dt + self.dl)
self.voltage[0] = h * (self.voltage[1] - self.abcv0) + self.abcv1 self.voltage[0] = h * (self.voltage[1] - self.abcv0) + self.abcv1
self.abcv0 = self.voltage[0] self.abcv0 = self.voltage[0]
self.abcv1 = self.voltage[1] self.abcv1 = self.voltage[1]
def update_voltage(self, time, G): def update_voltage(self, time, G):
"""Updates voltage values along the transmission line. """Updates voltage values along the transmission line.
Args: Args:
time (float): Absolute time. time (float): Absolute time.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
# Update all the voltage values along the line # Update all the voltage values along the line
self.voltage[1:self.nl] -= self.resistance * (c * G.dt / self.dl) * (self.current[1:self.nl] - self.current[0:self.nl - 1]) 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 # 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) 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) * waveform.amp * waveform.calculate_value(time - 0.5 * G.dt, G.dt)
@@ -279,12 +302,13 @@ class TransmissionLine(Source):
def update_current(self, time, G): def update_current(self, time, G):
"""Updates current values along the transmission line. """Updates current values along the transmission line.
Args: Args:
time (float): Absolute time. time (float): Absolute time.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
# Update all the current values along the line # Update all the current values along the line
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]) 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])
@@ -293,23 +317,25 @@ class TransmissionLine(Source):
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] += (c * G.dt / self.dl) * waveform.amp * waveform.calculate_value(time - 0.5 * G.dt, G.dt) * (1 / self.resistance)
def update_electric(self, abstime, Ex, Ey, Ez, G): def update_electric(self, abstime, Ex, Ey, Ez, G):
"""Updates electric field value in the main grid from voltage value in the transmission line. """Updates electric field value in the main grid from voltage value
in the transmission line.
Args: Args:
abstime (float): Absolute time. abstime (float): Absolute time.
Ex, Ey, Ez (memory view): numpy array of electric field values. Ex, Ey, Ez (memory view): numpy array of electric field values.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
if abstime >= self.start and abstime <= self.stop: if abstime >= self.start and abstime <= self.stop:
# Set the time of the waveform evaluation to account for any delay in the start # Set the time of the waveform evaluation to account
# for any delay in the start
time = abstime - self.start time = abstime - self.start
i = self.xcoord i = self.xcoord
j = self.ycoord j = self.ycoord
k = self.zcoord k = self.zcoord
self.update_voltage(time, G) self.update_voltage(time, G)
if self.polarisation is 'x': if self.polarisation is 'x':
Ex[i, j, k] = - self.voltage[self.antpos] / G.dx Ex[i, j, k] = - self.voltage[self.antpos] / G.dx
@@ -320,21 +346,24 @@ class TransmissionLine(Source):
Ez[i, j, k] = - self.voltage[self.antpos] / G.dz Ez[i, j, k] = - self.voltage[self.antpos] / G.dz
def update_magnetic(self, abstime, Hx, Hy, Hz, G): def update_magnetic(self, abstime, Hx, Hy, Hz, G):
"""Updates current value in transmission line from magnetic field values in the main grid. """Updates current value in transmission line from magnetic field
values in the main grid.
Args: Args:
abstime (float): Absolute time. abstime (float): Absolute time.
Hx, Hy, Hz (memory view): numpy array of magnetic field values. Hx, Hy, Hz (memory view): numpy array of magnetic field values.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters
describing the model.
""" """
if abstime >= self.start and abstime <= self.stop: if abstime >= self.start and abstime <= self.stop:
# Set the time of the waveform evaluation to account for any delay in the start # Set the time of the waveform evaluation to account
# for any delay in the start
time = abstime - self.start time = abstime - self.start
i = self.xcoord i = self.xcoord
j = self.ycoord j = self.ycoord
k = self.zcoord k = self.zcoord
if self.polarisation is 'x': if self.polarisation is 'x':
self.current[self.antpos] = Ix(i, j, k, G.Hy, G.Hz, G) self.current[self.antpos] = Ix(i, j, k, G.Hy, G.Hz, G)
@@ -346,3 +375,15 @@ class TransmissionLine(Source):
self.update_current(time, G) self.update_current(time, G)
class TEMTransmissionLine(TransmissionLine):
def __init__(self, G):
super().__init__(G)
self.type_str = 'TEMTransmissionLine'
def update_magnetic(self, abstime, Hx, Hy, Hz, G):
pass
def update_electric(self, abstime, Ex, Ey, Ez, G):
pass