publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-07 15:10:13 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Merge pull request #373 from AdittyaPal/devel

浏览代码 
[GSoC 2023] Implementing a discrete plane wave source
这个提交包含在：
Craig Warren
2025-03-10 13:38:49 +00:00
提交者 GitHub
父节点 eaef3c2127 4fe9fd52d7
当前提交 ac24757ebc
共有 10 个文件被更改，包括 1510 次插入 和 3 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

2
conda_env.yml
查看文件

@@ -1,7 +1,7 @@
# This file may be used to create a conda environment using:
# $ conda env create -f <this file>
name: gprMax-devel
name: gprMax-devel_PW
dependencies:
- python>3.7
- colorama

3
gprMax/__init__.py
查看文件

@@ -46,6 +46,8 @@ from .cmds_multiuse import (
VoltageSource,
Waveform,
)
from .cmds_multiuse import DiscretePlaneWave
from .cmds_singleuse import (
Discretisation,
Domain,
@@ -62,4 +64,5 @@ from .gprMax import run as run
from .scene import Scene
from .subgrids.user_objects import SubGridHSG
__name__ = "gprMax"

129
gprMax/cmds_multiuse.py
查看文件

@@ -18,7 +18,6 @@
import inspect
import logging
from pathlib import Path
import numpy as np
from scipy import interpolate
@@ -43,6 +42,7 @@ from .sources import HertzianDipole as HertzianDipoleUser
from .sources import MagneticDipole as MagneticDipoleUser
from .sources import TransmissionLine as TransmissionLineUser
from .sources import VoltageSource as VoltageSourceUser
from .sources import DiscretePlaneWave as DiscretePlaneWaveUser
from .subgrids.grid import SubGridBaseGrid
from .utilities.utilities import round_value
from .waveforms import Waveform as WaveformUser
@@ -946,6 +946,133 @@ class TransmissionLine(UserObjectMulti):
grid.transmissionlines.append(t)
"""
------------------------------------------------------------------------------
Add the UserMultiObject Class for the Discrete Plane Wave Implementation
------------------------------------------------------------------------------
"""
class DiscretePlaneWave(UserObjectMulti):
"""
Specifies a plane wave implemented using the discrete plane wave formulation.
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.order = 19
self.hash = "#discrete_plane_wave"
def create(self, grid, uip):
try:
theta = self.kwargs["theta"]
phi = self.kwargs["phi"]
psi = self.kwargs["psi"]
p1 = self.kwargs["p1"]
p2 = self.kwargs["p2"]
waveform_id = self.kwargs["waveform_id"]
except KeyError:
logger.exception(f"{self.params_str()} requires at least ten parameters.")
raise
try:
dtheta = self.kwargs["delta_theta"]
dphi = self.kwargs["delta_phi"]
except KeyError:
dtheta = 1.0
dphi = 1.0
# Warn about using a discrete plane wave on GPU
if config.sim_config.general["solver"] in ["cuda", "opencl"]:
logger.exception(
f"{self.params_str()} cannot currently be used "
+ "with the CUDA or OpenCL-based solver. "
)
raise ValueError
# Check if there is a waveformID in the waveforms list
if not any(x.ID == waveform_id for x in grid.waveforms):
logger.exception(f"{self.params_str()} there is no waveform " + f"with the identifier {waveform_id}.")
raise ValueError
if (theta > 180):
theta -= (np.floor(theta/180) * 180.0)
if (phi > 360):
phi -= (np.floor(phi/360) * 360.0)
if (psi > 360):
psi -= (np.floor(psi/360) * 360.0)
x_start, y_start, z_start = uip.check_src_rx_point(p1, self.params_str())
x_stop, y_stop, z_stop = uip.check_src_rx_point(p2, self.params_str())
DPW = DiscretePlaneWaveUser(grid)
DPW.corners = np.array([x_start, y_start, z_start, x_stop, y_stop, z_stop], dtype=np.int32)
DPW.waveformID = waveform_id
DPW.initializeDiscretePlaneWave(psi, phi, dphi, theta, dtheta, grid)
try:
DPW.material_ID = self.kwargs["material_id"]
except KeyError:
DPW.material_ID = 1
try:
# Check source start & source remove time parameters
start = self.kwargs["start"]
stop = self.kwargs["stop"]
if start < 0:
logger.exception(
self.params_str() + (" delay of the initiation " "of the source should not " "be less than zero."))
raise ValueError
if stop < 0:
logger.exception(self.params_str() + (" time to remove the " "source should not be " "less than zero."))
raise ValueError
if stop - start <= 0:
logger.exception(self.params_str() + (" duration of the source " "should not be zero or " "less."))
raise ValueError
DPW.start = start
DPW.stop = min(stop, grid.timewindow)
startstop = f" start time {t.start:g} secs, finish time " + f"{t.stop:g} secs "
except KeyError:
DPW.start = 0
DPW.stop = grid.timewindow
startstop = " "
precompute = True
if(precompute == True):
DPW.calculate_waveform_values(grid)
logger.info(
f"{self.grid_name(grid)} Discrete Plane Wave within the TFSF Box "
+ f"spanning from {p1} m to {p2} m, incident in the direction "
+ f"theta {theta} degrees and phi {phi} degrees "
+ startstop
+ f"using waveform {DPW.waveformID} created."
)
phi_approx = np.arctan2(DPW.m[1]/grid.dy, DPW.m[0]/grid.dx)*180/np.pi
theta_approx = np.arctan2( np.sqrt( (DPW.m[0]/grid.dx)*(DPW.m[0]/grid.dx) + (DPW.m[1]/grid.dy)*(DPW.m[1]/grid.dy) ),
DPW.m[2]/grid.dz) * 180/np.pi
logger.info(
f"Since {self.grid_name(grid)} Discrete Plane Wave has been discretized "
+ f"the angles have been approzimated to the nearest rational angles "
+ f"with some small tolerance levels. The chosen rational integers are "
+ f"[m_x, m_y, m_z] : {DPW.m[:3]}. The approximated angles are : \n"
+ "Approximated Phi : {:.3f}".format(phi_approx) + "\nApproximated Theta : {:.3f} .".format(theta_approx)
)
grid.discreteplanewaves.append(DPW)
'''
DPW = []
start = time.time()
for i in range(self.noOfWaves):
DPW.append(DiscretePlaneWaveUser(self.time_duration, self.dimensions, number_x, number_y, number_z))
SpaceGrid = TFSFBoxUser(number_x, number_y, number_z, self.corners, self.time_duration, self.dimensions, self.noOfWaves)
SpaceGrid.getFields(DPW, self.snapshot, angles, number, self.dx, self.dy, self.dz, self.dt, self.ppw)
'''
"""
------------------------------------------------------------------------------
End of the UserMultiObject Class for the Discrete Plane Wave Implementation
------------------------------------------------------------------------------
"""
class Rx(UserObjectMulti):
"""Specifies output points in the model.

904
gprMax/cython/planeWaveModules.pyx 普通文件
查看文件

@@ -0,0 +1,904 @@
import numpy as np
cimport cython
from libc.math cimport floor, ceil, round, sqrt, tan, cos, sin, atan2, abs, pow, exp, M_PI
from libc.stdio cimport FILE, fopen, fwrite, fclose
from libc.string cimport strcmp
from cython.parallel import prange
from gprMax.config cimport float_or_double
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cpdef double[:, ::1] getProjections(
double psi,
int[:] m
):
'''
Method to get the projection vectors to source the magnetic fields of the plane wave.
__________________________
Input parameters:
--------------------------
psi, float : the angle describing the polatan value of the required phi angle (which would be approximated to a rational number)
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles, for assignment of the correct element
to the three dimensional FDTD grid from the one dimensional representation, last element stores max(m_x, m_y, m_z)
__________________________
Returns:
--------------------------
projections, float array : stores the projections for sourcing the magnetic field and the sourcing vector
'''
cdef double phi, theta, cos_phi, sin_phi, cos_psi, sin_psi, cos_theta, sin_theta
cdef double[:, ::1] projections = np.zeros((2, 3), order='C')
if m[0] == 0:
phi = M_PI / 2
else:
phi = atan2(m[1], m[0])
if m[2] == 0:
theta = M_PI / 2
else:
theta = atan2(sqrt(m[0] * m[0] + m[1] * m[1]), m[2])
cos_phi = cos(phi)
sin_phi = sin(phi)
cos_theta = cos(theta)
sin_theta = sin(theta)
cos_psi = cos(psi)
sin_psi = sin(psi)
# Magnetic field projection vector
projections[0, 0] = sin_psi * sin_phi + cos_psi * cos_theta * cos_phi
projections[0, 1] = -sin_psi * cos_phi + cos_psi * cos_theta * sin_phi
projections[0, 2] = -cos_psi * sin_theta
# Direction cosines
projections[1, 0] = sin_theta * cos_phi
projections[1, 1] = sin_theta * sin_phi
projections[1, 2] = cos_theta
return projections
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cdef int[:] getPhi(
int[:, :] integers,
double required_ratio,
double tolerance
):
'''
Method to get the rational angle approximation to phi within the requested tolerance level using
Farey Fractions to determine a rational number closest to the real number.
__________________________
Input parameters:
--------------------------
integers, int array : array of integers to determine the value of m_x and m_y.
required_ratio, float : tan value of the required phi angle (which would be approximated to a rational number)
tolerance, float : acceptable deviation in the tan value of the rational angle from phi
__________________________
Returns:
--------------------------
integers, int array : sequence of the two integers [m_x, m_y]
'''
if(abs(integers[2, 0]/<double>integers[2, 1]-required_ratio)<tolerance):
return integers[2, :]
while True:
integers[1, 0] = integers[0, 0]+integers[2, 0]
integers[1, 1] = integers[0, 1]+integers[2, 1]
ratio = integers[1, 0]/<double>integers[1, 1]
if (abs(ratio - required_ratio) <= tolerance):
return integers[1, :]
elif(ratio >= required_ratio):
integers[2, :] = integers[1, :]
else:
integers[0, :] = integers[1, :]
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cdef inline double getTanValue(
int[:] integers,
double[:] dr
):
'''
Method to return the tan value of the angle approximated to theta given the three integers.
__________________________
Input parameters:
--------------------------
integers, int array : three integers for the rational angle approximation
dr, double array : array containing the separation between grid points along the three axes [dx, dy, dz]
__________________________
Returns:
--------------------------
_tanValue, double : tan value of the rationa angle cprrenponding to the integers m_x, m_y, m_z
'''
if(integers[2]==0): #if rational angle==90 degrees
return 99999.0 #return a large number to avoid division by zero error
else:
return sqrt((integers[0]/dr[0])*(integers[0]/dr[0]) + (integers[1]/dr[1])*(integers[1]/dr[1])
)/(integers[2]/dr[2])
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cdef int[:, :] get_mZ(
int m_x,
int m_y,
double theta,
double[:] Delta_r
):
'''
Method to get the arrays to perform a binary search to determine a rational number, m_z, closest to the real number, m_z, to get the desired tan Theta value.
__________________________
Input parameters:
--------------------------
m_x and m_y, int : integers approximating the rational angle to the tan value of phi
theta, float : polar angle of the incident plane wave (in radians) to be approximated to a rational angle
Delta_r, float array : projections of the propagation vector along the different coordinate axes
__________________________
Returns:
--------------------------
_integers, int array : a two dimensional sequence of the three integers [m_x, m_y, m_z] to perform a binary search
to determine the value of m_z within the given limits.
'''
cdef double m_z = 0
m_z = sqrt((m_x/Delta_r[0])*(m_x/Delta_r[0]) + (m_y/Delta_r[1])*(m_y/Delta_r[1]))/(tan(theta)/Delta_r[2]) #get an estimate of the m_z value
return np.array([[m_x, m_y, floor(m_z)],
[m_x, m_y, round(m_z)],
[m_x, m_y, ceil(m_z)]], dtype=np.int32, order='C') #set up the integer array to search for an appropriate m_z
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cdef int[:] getTheta(
int m_x,
int m_y,
double theta,
double Delta_theta,
double[:] Delta_r
):
'''
Method to get the rational angle approximation to theta within the requested tolerance level using
Binary Search to determine a rational number closest to the real number.
__________________________
Input parameters:
--------------------------
m_x and m_y, int : integers approximating the rational angle to the tan value of phi
theta, float : polar angle of the incident plane wave (in radians) to be approximated to a rational angle
Delta_theta, float : permissible error in the rational angle approximation to theta (in radians)
Delta_r, float array : projections of the propagation vector along the different coordinate axes
__________________________
Returns:
--------------------------
integers, int array : sequence of the three integers [m_x, m_y, m_z]
'''
cdef Py_ssize_t i, j = 0
cdef double tan_theta = 0.0
cdef int[:, :] integers = get_mZ(m_x, m_y, theta, Delta_r) #set up the integer array to search for an appropriate m_z
while True: #if tan value of m_z greater than permitted tolerance
tan_theta = getTanValue(integers[1, :], Delta_r)
if(abs(tan_theta - tan(theta)) <= Delta_theta / (cos(theta) * cos(theta))):
break
for i in range(3):
for j in range(3):
integers[i, j] = 2*integers[i, j] #expand the serach space by multiplying 2 to both the numerator and the denominator
while(integers[0, 2]<integers[2, 2]-1): #while there are integers to search for in the denominator
integers[1, 2] = (integers[0, 2]+integers[2, 2])/2 #get the mid integer between the upper and lower limits of denominator
tan_theta = getTanValue(integers[1, :], Delta_r)
if(tan_theta < tan(theta)): #if m_z results in a smaller tan value, make the denominator smaller
integers[2, 2] = integers[1, 2] #decrease m_z, serach in the lower half of the sample space
elif(tan_theta > tan(theta)): #if m_z results in a larger tan value, make the denominator larger
integers[0, 2] = integers[1, 2] #increase m_z, serach in the upper half of the sample space
return integers[1, :]
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cpdef int[:, ::1] getIntegerForAngles(
double phi,
double Delta_phi,
double theta,
double Delta_theta,
double[:] Delta_r
):
'''
Method to get [m_x, m_y, m_z] to determine the rational angles given phi and theta along with the permissible tolerances
__________________________
Input parameters:
--------------------------
phi, float : azimuthal angle of the incident plane wave (in degrees) to be approximated to a rational angle
Delta_phi, float : permissible error in the rational angle approximation to phi (in degrees)
theta, float : polar angle of the incident plane wave (in degrees) to be approximated to a rational angle
Delta_theta, float : permissible error in the rational angle approximation to theta (in degrees)
Delta_r, float : projections of the propagation vector along the different coordinate axes
__________________________
Returns:
--------------------------
quadrants[0, :], int array : the integers specifying the direction of propagation of the plane wave along the three coordinate axes
quadrants[1, :], int array : sequence of the three integers [m_x, m_y, m_z]
'''
cdef double required_ratio_phi, tolerance_phi = 0.0
cdef int m_x, m_y, m_z = 0
cdef int[:, ::1] quadrants = np.ones((2, 3), dtype=np.int32)
if(theta>=90):
quadrants[0, 2] = -1
theta = 180-theta
if(phi>=90 and phi<180):
quadrants[0, 0] = -1
phi = 180-phi
elif(phi>=180 and phi<270):
quadrants[0, 0] = -1
quadrants[0, 1] = -1
phi = phi-180
elif(phi>=270 and phi<360):
quadrants[0, 1] = -1
phi = 360-phi
if(0 <= phi < 90): #handle the case of phi=90 degrees separately
required_ratio_phi = tan(M_PI/180*phi) * Delta_r[1] / Delta_r[0] #to avoid division by zero exception
tolerance_phi = M_PI/180*Delta_phi / (cos(M_PI/180*phi)*cos(M_PI/180*phi)) * Delta_r[1] / Delta_r[0] #get the persissible error in tan phi
m_y, m_x = getPhi(np.array([[floor(required_ratio_phi), 1],
[required_ratio_phi, 1],
[ceil(required_ratio_phi), 1]], dtype=np.int32, order='C')
, required_ratio_phi, tolerance_phi) #get the integers [m_x, m_y] for rational angle approximation to phi
else:
m_x = 0
m_y = 1
if(theta < 90):
m_x, m_y, m_z = getTheta(m_x, m_y, M_PI/180*theta, M_PI/180*Delta_theta, Delta_r) #get the integers [m_x, m_y, m_z] for rational angle approximation to theta
else: #handle the case of theta=90 degrees separately
m_z = 0 #to avoid division by zero exception
quadrants[1, 0] = m_x
quadrants[1, 1] = m_y
quadrants[1, 2] = m_z
return quadrants
@cython.wraparound(False)
@cython.boundscheck(False)
cdef void applyTFSFMagnetic(
int nthreads,
float_or_double[:, :, ::1] Hx,
float_or_double[:, :, ::1] Hy,
float_or_double[:, :, ::1] Hz,
float_or_double[:, ::1] E_fields,
float_or_double[:] updatecoeffsH,
int[:] m,
int[:] corners
):
'''
Method to implement the total field-scattered field formulation for the magnetic field on the edge of the TF/SF region of the TFSF Box.
__________________________
Input parameters:
--------------------------
nthreads, int : number of threads to parallelize the for loops on
Hx, Hy, Hz, double array : stores the magnetic fields for the grid cells over the TFSF box at particular indices
E_fields, double array : array to store the electric fields of the one dimensional representation of the plane wave
in a direction along which the wave propagates
updatecoeffsH, float : stores the coefficients of the fields in the TFSF assignment equation for the magnetic field
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles,
for assignment of the correct element to the three dimensional FDTD grid
from the one dimensional representation, last element stores max(m_x, m_y, m_z)
corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries
'''
cdef Py_ssize_t i, j, k = 0
# Precompute index values
cdef int index = 0
cdef int m_x = m[0]
cdef int m_y = m[1]
cdef int m_z = m[2]
cdef int x_start = corners[0]
cdef int y_start = corners[1]
cdef int z_start = corners[2]
cdef int x_stop = corners[3]
cdef int y_stop = corners[4]
cdef int z_stop = corners[5]
cdef float_or_double[:] E_x = E_fields[0, :]
cdef float_or_double[:] E_y = E_fields[1, :]
cdef float_or_double[:] E_z = E_fields[2, :]
cdef float_or_double coef_H_xy = updatecoeffsH[2]
cdef float_or_double coef_H_xz = updatecoeffsH[3]
cdef float_or_double coef_H_yz = updatecoeffsH[3]
cdef float_or_double coef_H_yx = updatecoeffsH[1]
cdef float_or_double coef_H_zx = updatecoeffsH[1]
cdef float_or_double coef_H_zy = updatecoeffsH[2]
#**** constant x faces -- scattered-field nodes ****
i = x_start
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Hy at firstX-1/2 by subtracting Ez_inc
index = m_x * i + m_y * j + m_z * k
Hy[i-1, j, k] -= coef_H_yx * E_z[index]
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Hz at firstX-1/2 by adding Ey_inc
index = m_x * i + m_y * j + m_z * k
Hz[i-1, j, k] += coef_H_zx * E_y[index]
i = x_stop
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Hy at lastX+1/2 by adding Ez_inc
index = m_x * i + m_y * j + m_z * k
Hy[i, j, k] += coef_H_yx * E_z[index]
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Hz at lastX+1/2 by subtractinging Ey_inc
index = m_x * i + m_y * j + m_z * k
Hz[i, j, k] -= coef_H_zx * E_y[index]
#**** constant y faces -- scattered-field nodes ****
j = y_start
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Hx at firstY-1/2 by adding Ez_inc
index = m_x * i + m_y * j + m_z * k
Hx[i, j-1, k] += coef_H_xy * E_z[index]
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Hz at firstY-1/2 by subtracting Ex_inc
index = m_x * i + m_y * j + m_z * k
Hz[i, j-1, k] -= coef_H_zy * E_x[index]
j = y_stop
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Hx at lastY+1/2 by subtracting Ez_inc
index = m_x * i + m_y * j + m_z * k
Hx[i, j, k] -= coef_H_xy * E_z[index]
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Hz at lastY-1/2 by adding Ex_inc
index = m_x * i + m_y * j + m_z * k
Hz[i, j, k] += coef_H_zy * E_x[index]
#**** constant z faces -- scattered-field nodes ****
k = z_start
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop+1):
#correct Hy at firstZ-1/2 by adding Ex_inc
index = m_x * i + m_y * j + m_z * k
Hy[i, j, k-1] += coef_H_yz * E_x[index]
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop):
#correct Hx at firstZ-1/2 by subtracting Ey_inc
index = m_x * i + m_y * j + m_z * k
Hx[i, j, k-1] -= coef_H_xz * E_y[index]
k = z_stop
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop+1):
#correct Hy at firstZ-1/2 by subtracting Ex_inc
index = m_x * i + m_y * j + m_z * k
Hy[i, j, k] -= coef_H_yz * E_x[index]
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop):
#correct Hx at lastZ+1/2 by adding Ey_inc
index = m_x * i + m_y * j + m_z * k
Hx[i, j, k] += coef_H_xz * E_y[index]
@cython.wraparound(False)
@cython.boundscheck(False)
cdef void applyTFSFElectric(
int nthreads,
float_or_double[:, :, ::1] Ex,
float_or_double[:, :, ::1] Ey,
float_or_double[:, :, ::1] Ez,
float_or_double[:, ::1] H_fields,
float_or_double[:] updatecoeffsE,
int[:] m,
int[:] corners
):
'''
Method to implement the total field-scattered field formulation for the electric field on the edge of the TF/SF region of the TFSF Box.
__________________________
Input parameters:
--------------------------
nthreads, int : number of threads to parallelize the for loops on
Ex, Ey, Ez, double array : stores the magnetic fields for the grid cells over the TFSF box at particular indices
H_fields, double array : array to store the electric fields of the one dimensional representation of the plane wave
in a direction along which the wave propagates
updatecoeffsE, float : stores the coefficients of the fields in the TFSF assignment equation for the magnetic field
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles,
for assignment of the correct element to the three dimensional FDTD grid
from the one dimensional representation, last element stores max(m_x, m_y, m_z)
corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries
'''
cdef Py_ssize_t i, j, k = 0
# Precompute index values
cdef int index = 0
cdef int m_x = m[0]
cdef int m_y = m[1]
cdef int m_z = m[2]
cdef int x_start = corners[0]
cdef int y_start = corners[1]
cdef int z_start = corners[2]
cdef int x_stop = corners[3]
cdef int y_stop = corners[4]
cdef int z_stop = corners[5]
cdef float_or_double[:] H_x = H_fields[0, :]
cdef float_or_double[:] H_y = H_fields[1, :]
cdef float_or_double[:] H_z = H_fields[2, :]
cdef float_or_double coef_E_xy = updatecoeffsE[2]
cdef float_or_double coef_E_xz = updatecoeffsE[3]
cdef float_or_double coef_E_yz = updatecoeffsE[3]
cdef float_or_double coef_E_yx = updatecoeffsE[1]
cdef float_or_double coef_E_zx = updatecoeffsE[1]
cdef float_or_double coef_E_zy = updatecoeffsE[2]
#**** constant x faces -- total-field nodes ****/
i = x_start
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Ez at firstX face by subtracting Hy_inc
index = m_x * (i-1) + m_y * j + m_z * k
Ez[i, j, k] -= coef_E_zx * H_y[index]
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Ey at firstX face by adding Hz_inc
index = m_x * (i-1) + m_y * j + m_z * k
Ey[i, j, k] += coef_E_yx * H_z[index]
i = x_stop
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Ez at lastX face by adding Hy_inc
index = m_x * i + m_y * j + m_z * k
Ez[i, j, k] += coef_E_zx * H_y[index]
i = x_stop
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Ey at lastX face by subtracting Hz_inc
index = m_x * i + m_y * j + m_z * k
Ey[i, j, k] -= coef_E_yx * H_z[index]
#**** constant y faces -- total-field nodes ****/
j = y_start
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Ez at firstY face by adding Hx_inc
index = m_x * i + m_y * (j-1) + m_z * k
Ez[i, j, k] += coef_E_zy * H_x[index]
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Ex at firstY face by subtracting Hz_inc
index = m_x * i + m_y * (j-1) + m_z * k
Ex[i, j, k] -= coef_E_xy * H_z[index]
j = y_stop
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop):
#correct Ez at lastY face by subtracting Hx_inc
index = m_x * i + m_y * j + m_z * k
Ez[i, j, k] -= coef_E_zy * H_x[index]
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for k in range(z_start, z_stop+1):
#correct Ex at lastY face by adding Hz_inc
index = m_x * i + m_y * j + m_z * k
Ex[i, j, k] += coef_E_xy * H_z[index]
#**** constant z faces -- total-field nodes ****/
k = z_start
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop):
#correct Ey at firstZ face by subtracting Hx_inc
index = m_x * i + m_y * j + m_z * (k-1)
Ey[i, j, k] -= coef_E_yz * H_x[index]
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop+1):
#correct Ex at firstZ face by adding Hy_inc
index = m_x * i + m_y * j + m_z * (k-1)
Ex[i, j, k] += coef_E_xz * H_y[index]
k = z_stop
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop):
#correct Ey at lastZ face by adding Hx_inc
index = m_x * i + m_y * j + m_z * k
Ey[i, j, k] += coef_E_yz * H_x[index]
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
for j in range(y_start, y_stop+1):
#correct Ex at lastZ face by subtracting Hy_inc
index = m_x * i + m_y * j + m_z * k
Ex[i, j, k] -= coef_E_xz * H_y[index]
cdef void initializeMagneticFields(
int[:] m,
float_or_double[:, ::1] H_fields,
double[:] projections,
float_or_double[:, ::1] waveformvalues_wholedt,
bint precompute,
int iteration,
double dt,
double ds,
double c,
double start,
double stop,
double freq,
char* wavetype
):
'''
Method to initialize the first few grid points of the the source waveform
__________________________
Input parameters:
--------------------------
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles,
for assignment of the correct element to the three dimensional FDTD grid
from the one dimensional representation, last element stores max(m_x, m_y, m_z)
H_fields, double array : array to store the electric fields of the one dimensional representation of the plane wave
in a direction along which the wave propagates
projections, float array : stores the projections of the magnetoc fields along the different cartesian axes.
waveformvalues_wholedt, double array : stores the precomputed waveforms at each timestep to initialize the magnetic fields
precompute, boolean : stores whether the fields values have been precomputed or should be computed on the fly
iterations, int : stores the number of iterations in the simulation
dt, float : stores the timestep for the simulation
ds, float : stores the projection vector for sourcing the plane wave
c, float : stores the sped of light in the medium
start, float : stores the start time at which the source is placed in the TFSF grid
stop, float : stores the stop time at which the source is removed from the TFSF grid
freq, float : the frequency if the introduced wave which determines the grid points per wavelength for the wave source
wavetype, string : stores the type of waveform whose magnitude should be returned
'''
cdef Py_ssize_t r = 0
cdef double time_x, time_y, time_z = 0.0
if (precompute == True):
for r in range(m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
H_fields[0, r] = projections[0] * waveformvalues_wholedt[0, r]
H_fields[1, r] = projections[1] * waveformvalues_wholedt[1, r]
H_fields[2, r] = projections[2] * waveformvalues_wholedt[2, r]
else:
for r in range(m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
time_x = dt * iteration - (r+ (m[1]+m[2])*0.5) * ds/c
time_y = dt * iteration - (r+ (m[2]+m[0])*0.5) * ds/c
time_z = dt * iteration - (r+ (m[0]+m[1])*0.5) * ds/c
if (dt * iteration >= start and dt * iteration <= stop):
# Set the time of the waveform evaluation to account for any delay in the start
H_fields[0, r] = projections[0] * getSource(time_x-start, freq, wavetype, dt)
H_fields[1, r] = projections[1] * getSource(time_y-start, freq, wavetype, dt)
H_fields[2, r] = projections[2] * getSource(time_z-start, freq, wavetype, dt)
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cdef void updateMagneticFields(
int n,
float_or_double[:, ::1] H_fields,
float_or_double[:, ::1] E_fields,
float_or_double[:] updatecoeffsH,
int[:] m
):
'''
Method to update the magnetic fields for the next time step using
Equation 8 of DOI: 10.1109/LAWP.2009.2016851
__________________________
Input parameters:
--------------------------
n , int : stores the spatial length of the DPW array to update each length grid cell
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
E_fields, double array : stores the electric fields of the DPW till temporal index time
updatecoeffsH, double array : stores the coefficients of the fields in the update equation for the magnetic field
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles,
for assignment of the correct element to the three dimensional FDTD grid
from the one dimensional representation, last element stores max(m_x, m_y, m_z)
'''
cdef Py_ssize_t j = 0
cdef float_or_double[:] E_x = E_fields[0, :]
cdef float_or_double[:] E_y = E_fields[1, :]
cdef float_or_double[:] E_z = E_fields[2, :]
cdef float_or_double[:] H_x = H_fields[0, :]
cdef float_or_double[:] H_y = H_fields[1, :]
cdef float_or_double[:] H_z = H_fields[2, :]
cdef float_or_double coef_H_xt = updatecoeffsH[0]
cdef float_or_double coef_H_xy = updatecoeffsH[2]
cdef float_or_double coef_H_xz = updatecoeffsH[3]
cdef float_or_double coef_H_yt = updatecoeffsH[0]
cdef float_or_double coef_H_yz = updatecoeffsH[3]
cdef float_or_double coef_H_yx = updatecoeffsH[1]
cdef float_or_double coef_H_zt = updatecoeffsH[0]
cdef float_or_double coef_H_zx = updatecoeffsH[1]
cdef float_or_double coef_H_zy = updatecoeffsH[2]
cdef int m_x = m[0]
cdef int m_y = m[1]
cdef int m_z = m[2]
for j in range(m[3], n-m[3]): #loop to update the magnetic field at each spatial index
H_x[j] = coef_H_xt * H_x[j] + coef_H_xz * ( E_y[j+m_z] - E_y[j] ) - coef_H_xy * ( E_z[j+m_y] - E_z[j] ) #equation 8 of Tan, Potter paper
H_y[j] = coef_H_yt * H_y[j] + coef_H_yx * ( E_z[j+m_x] - E_z[j] ) - coef_H_yz * ( E_x[j+m_z] - E_x[j] ) #equation 8 of Tan, Potter paper
H_z[j] = coef_H_zt * H_z[j] + coef_H_zy * ( E_x[j+m_y] - E_x[j] ) - coef_H_zx * ( E_y[j+m_x] - E_y[j] ) #equation 8 of Tan, Potter paper
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cdef void updateElectricFields(
int n,
float_or_double[:, ::1] H_fields,
float_or_double[:, ::1] E_fields,
float_or_double[:] updatecoeffsE,
int[:] m
):
'''
Method to update the electric fields for the next time step using
Equation 9 of DOI: 10.1109/LAWP.2009.2016851
__________________________
Input parameters:
--------------------------
n , int : stores the spatial length of the DPW array to update each length grid cell
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
E_fields, double array : stores the electric fields of the DPW till temporal index time
updatecoeffsE, double array : stores the coefficients of the fields in the update equation for the electric field
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles,
for assignment of the correct element to the three dimensional FDTD grid
from the one dimensional representation, last element stores max(m_x, m_y, m_z)
'''
cdef Py_ssize_t j = 0
cdef float_or_double[:] E_x = E_fields[0, :]
cdef float_or_double[:] E_y = E_fields[1, :]
cdef float_or_double[:] E_z = E_fields[2, :]
cdef float_or_double[:] H_x = H_fields[0, :]
cdef float_or_double[:] H_y = H_fields[1, :]
cdef float_or_double[:] H_z = H_fields[2, :]
cdef float_or_double coef_E_xt = updatecoeffsE[0]
cdef float_or_double coef_E_xy = updatecoeffsE[2]
cdef float_or_double coef_E_xz = updatecoeffsE[3]
cdef float_or_double coef_E_yt = updatecoeffsE[0]
cdef float_or_double coef_E_yz = updatecoeffsE[3]
cdef float_or_double coef_E_yx = updatecoeffsE[1]
cdef float_or_double coef_E_zt = updatecoeffsE[0]
cdef float_or_double coef_E_zx = updatecoeffsE[1]
cdef float_or_double coef_E_zy = updatecoeffsE[2]
cdef int m_x = m[0]
cdef int m_y = m[1]
cdef int m_z = m[2]
for j in range(m[3], n-m[3]): #loop to update the electric field at each spatial index
E_x[j] = coef_E_xt * E_x[j] + coef_E_xz * ( H_z[j] - H_z[j-m_y] ) - coef_E_xy * ( H_y[j] - H_y[j-m_z] ) #equation 9 of Tan, Potter paper
E_y[j] = coef_E_yt * E_y[j] + coef_E_yx * ( H_x[j] - H_x[j-m_z] ) - coef_E_yz * ( H_z[j] - H_z[j-m_x] ) #equation 9 of Tan, Potter paper
E_z[j] = coef_E_zt * E_z[j] + coef_E_zy * ( H_y[j] - H_y[j-m_x] ) - coef_E_zx * ( H_x[j] - H_x[j-m_y] ) #equation 9 of Tan, Potter paper
@cython.cdivision(True)
cpdef double getSource(
double time,
double freq,
char* wavetype,
double dt
):
'''
Method to get the magnitude of the source field in the direction perpendicular to the propagation of the plane wave
__________________________
Input parameters:
--------------------------
time, float : time at which the magnitude of the source is calculated
freq, float : the frequency if the introduced wave which determines the grid points per wavelength for the wave source
wavetype, string : stores the type of waveform whose magnitude should be returned
dt, double : the time upto which the wave should exist in a impulse delta pulse
__________________________
Returns:
--------------------------
sourceMagnitude, double : magnitude of the source for the requested indices at the current time
'''
# Waveforms
if (strcmp(wavetype, "gaussian") == 0):
return exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
elif (strcmp(wavetype, "gaussiandot") == 0 or strcmp(wavetype, "gaussianprime") == 0):
return -4.0 * M_PI * M_PI * freq * (time * freq - 1.0
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
elif (strcmp(wavetype, "gaussiandotnorm") == 0):
return -2.0 * M_PI * (time * freq - 1.0
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0))) * exp(0.5)
elif (strcmp(wavetype, "gaussiandotdot") == 0 or strcmp(wavetype, "gaussiandoubleprime") == 0):
return (2.0 * M_PI * freq) * (2.0 * M_PI * freq) * (2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)) - 1.0
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
elif (strcmp(wavetype, "gaussiandotdotnorm") == 0):
return (2.0 * (M_PI *(time * freq - 1.0)) * (M_PI * (time * freq - 1.0)) - 1.0
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
elif (strcmp(wavetype, "ricker") == 0):
return (1.0 - 2.0 * (M_PI *(time * freq - 1.0)) * (M_PI * (time * freq - 1.0))
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0))) # define a Ricker wave source
elif (strcmp(wavetype, "sine") == 0):
if (time * freq <= 1):
return sin(2.0 * M_PI * freq * time)
else:
return 0.0
elif (strcmp(wavetype, "contsine") == 0):
return min(0.25 * time* freq, 1) * sin(2 * M_PI * time* freq)
elif (strcmp(wavetype, "impulse") == 0):
if (time < dt): # time < dt condition required to do impulsive magnetic dipole
return 1.0
else:
return 0.0
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
cpdef void calculate1DWaveformValues(
float_or_double[:, :, ::1] waveformvalues_wholedt,
int iterations,
int[:] m,
double dt,
double ds,
double c,
double start,
double stop,
double freq,
char* wavetype
):
'''
Method to precompute the source waveform values so that the initialization if faster, if requested
__________________________
Input parameters:
--------------------------
waveformvalues_wholedt, double array : stores the precomputed waveforms at each timestep to initialize the magnetic fields
iterations, int : stores the number of iterations in the simulation
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles,
for assignment of the correct element to the three dimensional FDTD grid
from the one dimensional representation, last element stores max(m_x, m_y, m_z)
dt, float : stores the timestep for the simulation
ds, float : stores the projection vector for sourcing the plane wave
c, float : stores the sped of light in the medium
start, float : stores the start time at which the source is placed in the TFSF grid
stop, float : stores the stop time at which the source is removed from the TFSF grid
freq, float : the frequency if the introduced wave which determines the grid points per wavelength for the wave source
wavetype, string : stores the type of waveform whose magnitude should be returned
'''
cdef double time_x, time_y, time_z = 0.0
cdef Py_ssize_t iteration, r = 0
for iteration in range(iterations):
for r in range(m[3]):
time_x = dt * iteration - (r+ (m[1]+m[2])*0.5) * ds/c
time_y = dt * iteration - (r+ (m[2]+m[0])*0.5) * ds/c
time_z = dt * iteration - (r+ (m[0]+m[1])*0.5) * ds/c
if (dt * iteration >= start and dt * iteration <= stop):
# Set the time of the waveform evaluation to account for any delay in the start
waveformvalues_wholedt[iteration, 0, r] = getSource(time_x-start, freq, wavetype, dt)
waveformvalues_wholedt[iteration, 1, r] = getSource(time_y-start, freq, wavetype, dt)
waveformvalues_wholedt[iteration, 2, r] = getSource(time_z-start, freq, wavetype, dt)
cpdef void updatePlaneWave(
int n,
int nthreads,
float_or_double[:, ::1] H_fields,
float_or_double[:, ::1] E_fields,
float_or_double[:] updatecoeffsE,
float_or_double[:] updatecoeffsH,
float_or_double[:, :, ::1] Ex,
float_or_double[:, :, ::1] Ey,
float_or_double[:, :, ::1] Ez,
float_or_double[:, :, ::1] Hx,
float_or_double[:, :, ::1] Hy,
float_or_double[:, :, ::1] Hz,
double[:] projections,
float_or_double[:, ::1] waveformvalues_wholedt,
int[:] m,
int[:] corners,
bint precompute,
int iteration,
double dt,
double ds,
double c,
double start,
double stop,
double freq,
char* wavetype
):
initializeMagneticFields(m, H_fields, projections, waveformvalues_wholedt, precompute, iteration, dt, ds, c, start, stop, freq, wavetype)
updateMagneticFields(n, H_fields, E_fields, updatecoeffsH, m)
applyTFSFMagnetic(nthreads, Hx, Hy, Hz, E_fields, updatecoeffsH, m, corners)
applyTFSFElectric(nthreads, Ex, Ey, Ez, H_fields, updatecoeffsE, m, corners)
updateElectricFields(n, H_fields, E_fields, updatecoeffsE, m)
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef void takeSnapshot3D(double[:, :, ::1] field, char* filename):
'''
Method to write the fields of the plane wave simulation at a particular time step.
__________________________
Input parameters:
--------------------------
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices of the TFSF box at a particular time step
filename, string : specifies the file location where the fields are to be written
'''
cdef FILE *fptr = fopen(filename, "wb")
fwrite(&field[0, 0, 0], sizeof(double), field.size, fptr)
fclose(fptr)

1
gprMax/grid.py
查看文件

@@ -75,6 +75,7 @@ class FDTDGrid:
self.hertziandipoles = []
self.magneticdipoles = []
self.transmissionlines = []
self.discreteplanewaves = []
self.rxs = []
self.srcsteps = [0, 0, 0]
self.rxsteps = [0, 0, 0]

1
gprMax/hash_cmds_file.py
查看文件

@@ -252,6 +252,7 @@ def check_cmd_names(processedlines, checkessential=True):
"#hertzian_dipole",
"#magnetic_dipole",
"#transmission_line",
"#discrete_plane_wave",
"#excitation_file",
"#rx",
"#rx_array",

86
gprMax/hash_cmds_multiuse.py
查看文件

@@ -36,6 +36,7 @@ from .cmds_multiuse import (
Snapshot,
SoilPeplinski,
TransmissionLine,
DiscretePlaneWave,
VoltageSource,
Waveform,
)
@@ -220,6 +221,89 @@ def process_multicmds(multicmds):
scene_objects.append(tl)
cmdname = "#discrete_plane_wave"
if multicmds[cmdname] is not None:
for cmdinstance in multicmds[cmdname]:
tmp = cmdinstance.split()
if len(tmp) < 10:
logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least ten parameters")
raise ValueError
if len(tmp) == 10:
plWave = DiscretePlaneWave(
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
psi = float(tmp[6]),
phi = float(tmp[7]),
theta = float(tmp[8]),
waveform_id = tmp[9],
)
elif len(tmp) == 11:
plWave = DiscretePlaneWave(
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
psi = float(tmp[6]),
phi = float(tmp[7]),
theta = float(tmp[9]),
waveform_id = tmp[11],
material_ID = int(tmp[12]),
)
elif len(tmp) == 12:
plWave = DiscretePlaneWave(
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
psi = float(tmp[6]),
phi = float(tmp[7]),
delta_phi = float(tmp[8]),
theta = float(tmp[9]),
delta_theta = float(tmp[10]),
waveform_id = tmp[11],
)
elif len(tmp) == 13:
plWave = DiscretePlaneWave(
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
psi = float(tmp[6]),
phi = float(tmp[7]),
delta_phi = float(tmp[8]),
theta = float(tmp[9]),
delta_theta = float(tmp[10]),
waveform_id = tmp[11],
material_ID = int(tmp[12]),
)
elif len(tmp) == 14:
plWave = DiscretePlaneWave(
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
psi = float(tmp[6]),
phi = float(tmp[7]),
delta_phi = float(tmp[8]),
theta = float(tmp[9]),
delta_theta = float(tmp[10]),
waveform_id = tmp[11],
start = float(tmp[12]),
stop = float(tmp[13]),
)
elif len(tmp) == 15:
plWave = DiscretePlaneWave(
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
psi = float(tmp[6]),
phi = float(tmp[7]),
delta_phi = float(tmp[8]),
theta = float(tmp[9]),
delta_theta = float(tmp[10]),
waveform_id = tmp[11],
material_ID = int(tmp[12]),
start = float(tmp[13]),
stop = float(tmp[14]),
)
else:
logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters")
raise ValueError
scene_objects.append(plWave)
cmd = "#excitation_file"
if multicmds[cmd] is not None:
for cmdinstance in multicmds[cmdname]:
@@ -236,7 +320,7 @@ def process_multicmds(multicmds):
ex_file = ExcitationFile(filepath=tmp[0])
scene_objects.append(ex_file)
cmdname = "#rx"
if multicmds[cmdname] is not None:
for cmdinstance in multicmds[cmdname]:

1
gprMax/solvers.py
查看文件

@@ -112,6 +112,7 @@ class Solver:
self.updates.update_magnetic()
self.updates.update_magnetic_pml()
self.updates.update_magnetic_sources()
#self.updates.update_plane_waves()
if self.hsg:
self.updates.hsg_2()
self.updates.update_electric_a()

364
gprMax/sources.py
查看文件

@@ -25,6 +25,8 @@ import gprMax.config as config
from .fields_outputs import Ix, Iy, Iz
from .utilities.utilities import round_value
import scipy.constants as constants
from .cython.planeWaveModules import getIntegerForAngles, getProjections, calculate1DWaveformValues, updatePlaneWave, getSource
class Source:
"""Super-class which describes a generic source."""
@@ -642,3 +644,365 @@ class TransmissionLine(Source):
self.current[self.antpos] = Iz(i, j, k, G.Hx, G.Hy, G.Hz, G)
self.update_current(iteration, G)
class DiscretePlaneWave(Source):
"""
Class to implement the discrete plane wave (DPW) formulation as described in
Tan, T.; Potter, M. (2010).
FDTD Discrete Planewave (FDTD-DPW) Formulation for a Perfectly Matched Source in TFSF Simulations. ,
58(8), 0–2648. doi:10.1109/tap.2010.2050446
__________________________
Instance variables:
--------------------------
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles
last element stores max(m_x, m_y, m_z)
directions, int array : stores the directions of propagation of the DPW
dimensions, int : stores the number of dimensions in which the simulation is run (2D or 3D)
time_dimension, int : stores the time length over which the simulation is run
E_fileds, double array : stores the electric flieds associated with the 1D DPW
H_fields, double array : stores the magnetic fields associated with the 1D DPW
"""
def __init__(self, G):
"""
Defines the instance variables of class DiscretePlaneWave()
__________________________
Input parameters:
--------------------------
time_dimension, int : local variable to store the time length over which the simulation is run
dimensions, int : local variable to store the number of dimensions in which the simulation is run
"""
super().__init__()
self.m = np.zeros(3+1, dtype=np.int32) #+1 to store the max(m_x, m_y, m_z)
self.directions = np.zeros(3, dtype=np.int32)
self.length = 0
self.projections = np.zeros(3, dtype=config.sim_config.dtypes["float_or_double"])
self.corners = None
self.materialID = 1
self.ds = 0
def initializeDiscretePlaneWave(self, psi, phi, Delta_phi, theta, Delta_theta, G):
"""
Method to create a DPW, assign memory to the grids and get field values at different time and space indices
__________________________
Input parameters:
--------------------------
psi, float : polarization angle of the incident plane wave (in radians)
phi, float : azimuthal angle of the incident plane wave (in radians)
Delta_phi, float : permissible error in the rational angle approximation to phi (in radians)
theta, float : polar angle of the incident plane wave (in radians)
Delta_theta, float : permissible error in the rational angle approximation to theta (in radians)
number, int : number of cells in the 3D FDTD simulation
dx, double : separation between adjacent cells in the x direction
dy, double : separation between adjacent cells in the y direction
dz, double : separation between adjacent cells in the z direction
dt, double : time step for the FDTD simulation
__________________________
Returns:
--------------------------
E_fields, double array : the electric field for the DPW as it evolves over time and space indices
H_fields, double array : the magnetic field for the DPW as it evolves over time and space indices
C, double array : stores the coefficients of the fields for the update equation of the electric fields
D, double array : stores the coefficients of the fields for the update equation of the magnetic fields
"""
self.directions, self.m[:3] = getIntegerForAngles(phi, Delta_phi, theta, Delta_theta,
np.array([G.dx, G.dy, G.dz])) #get the integers for the nearest rational angle
self.m[3] = max(self.m) #store max(m_x, m_y, m_z) in the last element of the array
self.length = int( 2 * max(self.m[:-1]) * G.iterations) #set an appropriate length fo the one dimensional arrays
self.E_fields = np.zeros((3, self.length), order='C', dtype=config.sim_config.dtypes["float_or_double"])
self.H_fields = np.zeros((3, self.length), order='C', dtype=config.sim_config.dtypes["float_or_double"])
#the 1D grid has no ABC to terminate it, sufficiently long array prevents reflections from the back
# Projections for field components
projections_h, P = getProjections(psi*180/np.pi, self.m) #get the projection vertors for different fields
self.projections = projections_h / np.sqrt(config.m0/config.e0) #scale the projection vector for the mangetic field
if self.m[0] == 0: #calculate dr that is needed for sourcing the 1D array
if self.m[1] == 0:
if self.m[2] == 0:
raise ValueError("not all m_i values can be zero")
else:
self.ds = P[2]*G.dz/self.m[2]
else:
self.ds = P[1]*G.dy/self.m[1]
else:
self.ds = P[0]*G.dx/self.m[0]
def calculate_waveform_values(self, G, cythonize=True):
"""
Calculates all waveform values for source for duration of simulation.
Args:
G: FDTDGrid class describing a grid in a model.
"""
# Waveform values for sources that need to be calculated on whole timesteps
self.waveformvalues_wholedt = np.zeros((G.iterations, 3, self.m[3]), dtype=config.sim_config.dtypes["float_or_double"])
# Waveform values for sources that need to be calculated on half timesteps
#self.waveformvalues_halfdt = np.zeros((G.iterations), dtype=config.sim_config.dtypes["float_or_double"])
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
if(cythonize == True):
calculate1DWaveformValues(self.waveformvalues_wholedt, G.iterations, self.m, G.dt, self.ds, config.c,
self.start, self.stop, waveform.freq, waveform.type.encode('UTF-8'))
else:
for dimension in range(3):
for iteration in range(G.iterations):
for r in range(self.m[3]):
time = G.dt * iteration - (r+ (self.m[(dimension+1)%3]+self.m[(dimension+2)%3])*0.5) * self.ds/config.c
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.waveformvalues_wholedt[iteration, dimension, r] = waveform.calculate_value(time, G.dt)
#self.waveformvalues_halfdt[iteration] = waveform.calculate_value(time + 0.5 * G.dt, G.dt)
def update_plane_wave(self, nthreads, updatecoeffsE, updatecoeffsH, Ex, Ey, Ez, Hx, Hy, Hz, iteration, G, cythonize=True, precompute=True):
if(cythonize == True):
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
updatePlaneWave(self.length, nthreads, self.H_fields, self.E_fields,
updatecoeffsE[self.material_ID, :], updatecoeffsH[self.material_ID, :],
Ex, Ey, Ez, Hx, Hy, Hz, self.projections, self.waveformvalues_wholedt[iteration, :, :],
self.m, self.corners, precompute, iteration, G.dt, self.ds, config.c, self.start, self.stop,
waveform.freq, waveform.type.encode('UTF-8'))
else:
self.update_magnetic_field_1D(G, precompute)
self.apply_TFSF_conditions_magnetic(G)
self.apply_TFSF_conditions_electric(G)
self.update_electric_field_1D(G)
def initialize_magnetic_fields_1D(self, G, precompute):
if(precompute == True):
for dimension in range(3):
for r in range(self.m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
self.H_fields[dimension, r] = self.projections[dimension] * self.waveformvalues_wholedt[G.iteration, dimension, r]
#self.getSource(self.real_time - (j+(self.m[(i+1)%3]+self.m[(i+2)%3])*0.5)*self.ds/config.c)#, self.waveformID, G.dt)
else:
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
for dimension in range(3):
for r in range(self.m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
self.H_fields[dimension, r] = self.projections[dimension] * getSource(G.iteration*G.dt - (r+(self.m[(dimension+1)%3]+self.m[(dimension+2)%3])*0.5)*self.ds/config.c, waveform.freq, waveform.type.encode('UTF-8'), G.dt)
def update_magnetic_field_1D(self, G, precompute=True):
"""
Method to update the magnetic fields for the next time step using
Equation 8 of DOI: 10.1109/LAWP.2009.2016851
__________________________
Input parameters:
--------------------------
n , int : stores the spatial length of the DPW array so that each length grid cell is updated when the method updateMagneticFields() is called
H_coefficients, double array : stores the coefficients of the fields in the update equation for the magnetic field
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
E_fields, double array : stores the electric fields of the DPW till temporal index time
time, int : time index storing the current axis number which would be updated for the H_fields
__________________________
Returns:
--------------------------
H_fields, double array : magnetic field array with the axis entry for the current time added
"""
self.initialize_magnetic_fields_1D(G, precompute)
for i in range(3): #loop to update each component of the magnetic field
materialH = G.ID[3+i, (self.corners[0]+self.corners[3])//2,
(self.corners[1]+self.corners[4])//2,
(self.corners[2]+self.corners[5])//2 ]
for j in range(self.m[-1], self.length-self.m[-1]): #loop to update the magnetic field at each spatial index
self.H_fields[i, j] = G.updatecoeffsH[materialH, 0] * self.H_fields[i, j] + G.updatecoeffsH[materialH, (i+2)%3+1] * (
self.E_fields[(i+1)%3, j+self.m[(i+2)%3]] - self.E_fields[(i+1)%3, j]) - G.updatecoeffsH[materialH, (i+1)%3+1] * (
self.E_fields[(i+2)%3, j+self.m[(i+1)%3]] - self.E_fields[(i+2)%3, j]) #equation 8 of Tan, Potter paper
def update_electric_field_1D(self, G):
"""
Method to update the electric fields for the next time step using
Equation 9 of DOI: 10.1109/LAWP.2009.2016851
__________________________
Input parameters:
--------------------------
n , int : stores the spatial length of the DPW array so that each length grid cell is updated when the method updateMagneticFields() is called
E_coefficients, double array : stores the coefficients of the fields in the update equation for the electric field
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
E_fields, double array : stores the electric fields of the DPW till temporal index time
time, int : time index storing the current axis number which would be updated for the H_fields
__________________________
Returns:
--------------------------
E_fields, double array : electric field array with the axis entry for the current time added
"""
for i in range(3): #loop to update each component of the electric field
materialE = G.ID[i, (self.corners[0]+self.corners[3])//2,
(self.corners[1]+self.corners[4])//2,
(self.corners[2]+self.corners[5])//2 ]
for j in range(self.m[-1], self.length-self.m[-1]): #loop to update the electric field at each spatial index
self.E_fields[i, j] = G.updatecoeffsE[materialE, 0] * self.E_fields[i, j] + G.updatecoeffsE[materialE, (i+2)%3+1] * (
self.H_fields[(i+2)%3, j] - self.H_fields[(i+2)%3, j-self.m[(i+1)%3]]) - G.updatecoeffsE[materialE, (i+1)%3+1] * (
self.H_fields[(i+1)%3, j] - self.H_fields[(i+1)%3, j-self.m[(i+2)%3]]) #equation 9 of Tan, Potter paper
def getField(self, i, j, k, array, m, component):
return array[component, np.dot(m[:-1], np.array([i, j, k]))]
def apply_TFSF_conditions_magnetic(self, G):
#**** constant x faces -- scattered-field nodes ****
i = self.corners[0]
for j in range(self.corners[1], self.corners[4]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Hy at firstX-1/2 by subtracting Ez_inc
G.Hy[i-1, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 2)
for j in range(self.corners[1], self.corners[4]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Hz at firstX-1/2 by adding Ey_inc
G.Hz[i-1, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 1)
i = self.corners[3]
for j in range(self.corners[1], self.corners[4]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Hy at lastX+1/2 by adding Ez_inc
G.Hy[i, j, k] += G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 2)
for j in range(self.corners[1], self.corners[4]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Hz at lastX+1/2 by subtractinging Ey_inc
G.Hz[i, j, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 1)
#**** constant y faces -- scattered-field nodes ****
j = self.corners[1]
for i in range(self.corners[0], self.corners[3]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Hx at firstY-1/2 by adding Ez_inc
G.Hx[i, j-1, k] += G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 2)
for i in range(self.corners[0], self.corners[3]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Hz at firstY-1/2 by subtracting Ex_inc
G.Hz[i, j-1, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 0)
j = self.corners[4]
for i in range(self.corners[0], self.corners[3]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Hx at lastY+1/2 by subtracting Ez_inc
G.Hx[i, j, k] -= G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 2)
for i in range(self.corners[0], self.corners[3]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Hz at lastY-1/2 by adding Ex_inc
G.Hz[i, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 0)
#**** constant z faces -- scattered-field nodes ****
k = self.corners[2]
for i in range(self.corners[0], self.corners[3]):
for j in range(self.corners[1], self.corners[4]+1):
#correct Hy at firstZ-1/2 by adding Ex_inc
G.Hy[i, j, k-1] += G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 0)
for i in range(self.corners[0], self.corners[3]+1):
for j in range(self.corners[1], self.corners[4]):
#correct Hx at firstZ-1/2 by subtracting Ey_inc
G.Hx[i, j, k-1] -= G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 1)
k = self.corners[5]
for i in range(self.corners[0], self.corners[3]):
for j in range(self.corners[1], self.corners[4]+1):
#correct Hy at firstZ-1/2 by subtracting Ex_inc
G.Hy[i, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 0)
for i in range(self.corners[0], self.corners[3]+1):
for j in range(self.corners[1], self.corners[4]):
#correct Hx at lastZ+1/2 by adding Ey_inc
G.Hx[i, j, k] += G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 1)
def apply_TFSF_conditions_electric(self, G):
#**** constant x faces -- total-field nodes ****/
i = self.corners[0]
for j in range(self.corners[1], self.corners[4]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Ez at firstX face by subtracting Hy_inc
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(i-1, j, k, self.H_fields, self.m, 1)
for j in range(self.corners[1], self.corners[4]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Ey at firstX face by adding Hz_inc
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(i-1, j, k, self.H_fields, self.m, 2)
i = self.corners[3]
for j in range(self.corners[1], self.corners[4]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Ez at lastX face by adding Hy_inc
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(i, j, k, self.H_fields, self.m, 1)
i = self.corners[3]
for j in range(self.corners[1], self.corners[4]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Ey at lastX face by subtracting Hz_inc
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(i, j, k, self.H_fields, self.m, 2)
#**** constant y faces -- total-field nodes ****/
j = self.corners[1]
for i in range(self.corners[0], self.corners[3]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Ez at firstY face by adding Hx_inc
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(i, j-1, k, self.H_fields, self.m, 0)
for i in range(self.corners[0], self.corners[3]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Ex at firstY face by subtracting Hz_inc
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(i, j-1, k, self.H_fields, self.m, 2)
j = self.corners[4]
for i in range(self.corners[0], self.corners[3]+1):
for k in range(self.corners[2], self.corners[5]):
#correct Ez at lastY face by subtracting Hx_inc
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(i, j, k, self.H_fields, self.m, 0)
for i in range(self.corners[0], self.corners[3]):
for k in range(self.corners[2], self.corners[5]+1):
#correct Ex at lastY face by adding Hz_inc
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(i, j, k, self.H_fields, self.m, 2)
#**** constant z faces -- total-field nodes ****/
k = self.corners[2]
for i in range(self.corners[0], self.corners[3]+1):
for j in range(self.corners[1], self.corners[4]):
#correct Ey at firstZ face by subtracting Hx_inc
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(i, j, k-1, self.H_fields, self.m, 0)
for i in range(self.corners[0], self.corners[3]):
for j in range(self.corners[1], self.corners[4]+1):
#correct Ex at firstZ face by adding Hy_inc
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(i, j, k-1, self.H_fields, self.m, 1)
k = self.corners[5]
for i in range(self.corners[0], self.corners[3]+1):
for j in range(self.corners[1], self.corners[4]):
#correct Ey at lastZ face by adding Hx_inc
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(i, j, k, self.H_fields, self.m, 0)
for i in range(self.corners[0], self.corners[3]):
for j in range(self.corners[1], self.corners[4]+1):
#correct Ex at lastZ face by subtracting Hy_inc
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(i, j, k, self.H_fields, self.m, 1)

22
gprMax/updates.py
查看文件

@@ -101,6 +101,28 @@ class CPUUpdates:
self.grid.Hz,
self.grid,
)
"""Updates the magnetic and electric field components for the discrete plane wave."""
for source in self.grid.discreteplanewaves:
source.update_plane_wave(
config.get_model_config().ompthreads,
self.grid.updatecoeffsE,
self.grid.updatecoeffsH,
self.grid.Ex,
self.grid.Ey,
self.grid.Ez,
self.grid.Hx,
self.grid.Hy,
self.grid.Hz,
self.grid.iteration,
self.grid,
cythonize = True,
precompute = False,
)
if(self.grid.iteration % 10 == 0):
np.save('./snapshots/electric_z_{}.npy'.format(self.grid.iteration), self.grid.Ex)
def update_electric_a(self):
"""Updates electric field components."""