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

Updates to deal with Windows handling of complex values.

浏览代码
这个提交包含在：
Craig Warren
2022-01-18 14:39:08 +00:00
父节点 27b44ecc33
当前提交 8b1f419faf
共有 2 个文件被更改，包括 141 次插入 和 60 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

160
gprMax/templates/fields_updates_dispersive_template
查看文件

@@ -21,9 +21,16 @@ cimport numpy as np
from cython.parallel import prange
# Use C-functions from 'complex.h' but doesn't work for Windows as it doesn't
# support 'double complex' but instead defines its own type '_Dcomplex'.
# https://docs.microsoft.com/en-us/cpp/c-runtime-library/complex-math-support?view=vs-2019
# https://stackoverflow.com/questions/57837255/defining-dcomplex-externally-in-cython?rq=1
{% if not functions[0].iswin %}
cdef extern from "complex.h" nogil:
double creal(double complex z)
float crealf(float complex z)
{% endif %}
###############################################################
@@ -50,13 +57,15 @@ cpdef void {{ item.name_a }}(
{{ item.field_type }}[:, :, ::1] Hy,
{{ item.field_type }}[:, :, ::1] Hz
):
"""This function updates the electric field components when dispersive materials (with multiple poles) are present.
"""This function updates the electric field components when dispersive
materials (with multiple poles) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
maxpoles (int): Maximum number of poles
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
nx, ny, nz: int for grid size in cells.
nthreads: int for number of threads to use.
maxpoles: int for maximum number of poles.
updatecoeffs, T, ID, E, H: memoryviews to access to update coeffients,
temporary, ID and field component arrays.
"""
cdef Py_ssize_t i, j, k, pole
@@ -72,12 +81,22 @@ cpdef void {{ item.name_a }}(
phi = 0
for pole in range(maxpoles):
{% if 'complex' in item.dispersive_type %}
phi = phi + {{ item.real_part }}(updatecoeffsdispersive[material, pole * 3]) * {{ item.real_part }}(Tx[pole, i, j, k])
{% if item.iswin %}
phi = (phi + updatecoeffsdispersive[material, pole * 3].real
* Tx[pole, i, j, k].real)
{% else %}
phi = (phi + {{ item.real_part }}(updatecoeffsdispersive[material, pole * 3])
* {{ item.real_part }}(Tx[pole, i, j, k]))
{% endif %}
{% else %}
phi = phi + updatecoeffsdispersive[material, pole * 3] * Tx[pole, i, j, k]
{% endif %}
Tx[pole, i, j, k] = updatecoeffsdispersive[material, 1 + (pole * 3)] * Tx[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)] * Ex[i, j, k]
Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi
Tx[pole, i, j, k] = (updatecoeffsdispersive[material, 1 + (pole * 3)]
* Tx[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)]
* Ex[i, j, k])
Ex[i, j, k] = (updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2]
* (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3]
* (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi)
# Ey component
if nx != 1 or nz != 1:
@@ -88,12 +107,22 @@ cpdef void {{ item.name_a }}(
phi = 0
for pole in range(maxpoles):
{% if 'complex' in item.dispersive_type %}
phi = phi + {{ item.real_part }}(updatecoeffsdispersive[material, pole * 3]) * {{ item.real_part }}(Ty[pole, i, j, k])
{% if item.iswin %}
phi = (phi + updatecoeffsdispersive[material, pole * 3].real
* Ty[pole, i, j, k].real)
{% else %}
phi = (phi + {{ item.real_part }}(updatecoeffsdispersive[material, pole * 3])
* {{ item.real_part }}(Ty[pole, i, j, k]))
{% endif %}
{% else %}
phi = phi + updatecoeffsdispersive[material, pole * 3] * Ty[pole, i, j, k]
{% endif %}
Ty[pole, i, j, k] = updatecoeffsdispersive[material, 1 + (pole * 3)] * Ty[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)] * Ey[i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
Ty[pole, i, j, k] = (updatecoeffsdispersive[material, 1 + (pole * 3)]
* Ty[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)]
* Ey[i, j, k])
Ey[i, j, k] = (updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3]
* (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1]
* (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi)
# Ez component
if nx != 1 or ny != 1:
@@ -104,12 +133,22 @@ cpdef void {{ item.name_a }}(
phi = 0
for pole in range(maxpoles):
{% if 'complex' in item.dispersive_type %}
phi = phi + {{ item.real_part }}(updatecoeffsdispersive[material, pole * 3]) * {{ item.real_part }}(Tz[pole, i, j, k])
{% if item.iswin %}
phi = (phi + updatecoeffsdispersive[material, pole * 3].real
* Tz[pole, i, j, k].real)
{% else %}
phi = (phi + {{ item.real_part }}(updatecoeffsdispersive[material, pole * 3])
* {{ item.real_part }}(Tz[pole, i, j, k]))
{% endif %}
{% else %}
phi = phi + updatecoeffsdispersive[material, pole * 3] * Tz[pole, i, j, k]
{% endif %}
Tz[pole, i, j, k] = updatecoeffsdispersive[material, 1 + (pole * 3)] * Tz[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)] * Ez[i, j, k]
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi
Tz[pole, i, j, k] = (updatecoeffsdispersive[material, 1 + (pole * 3)]
* Tz[pole, i, j, k] + updatecoeffsdispersive[material, 2 + (pole * 3)]
* Ez[i, j, k])
Ez[i, j, k] = (updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1]
* (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2]
* (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi)
{% endfor %}
@@ -133,13 +172,15 @@ cpdef void {{ item.name_b }}(
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez
):
"""This function updates a temporary dispersive material array when disperisive materials (with multiple poles) are present.
"""This function updates a temporary dispersive material array when
disperisive materials (with multiple poles) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
maxpoles (int): Maximum number of poles
updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
nx, ny, nz: int for grid size in cells.
nthreads: int for number of threads to use.
maxpoles: int for maximum number of poles.
updatecoeffs, T, ID, E, H: memoryviews to access to update coeffients,
temporary, ID and field component arrays.
"""
cdef Py_ssize_t i, j, k, pole
@@ -152,7 +193,9 @@ cpdef void {{ item.name_b }}(
for k in range(1, nz):
material = ID[0, i, j, k]
for pole in range(maxpoles):
Tx[pole, i, j, k] = Tx[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ex[i, j, k]
Tx[pole, i, j, k] = (Tx[pole, i, j, k]
- updatecoeffsdispersive[material, 2 + (pole * 3)]
* Ex[i, j, k])
# Ey component
if nx != 1 or nz != 1:
@@ -161,7 +204,9 @@ cpdef void {{ item.name_b }}(
for k in range(1, nz):
material = ID[1, i, j, k]
for pole in range(maxpoles):
Ty[pole, i, j, k] = Ty[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ey[i, j, k]
Ty[pole, i, j, k] = (Ty[pole, i, j, k]
- updatecoeffsdispersive[material, 2 + (pole * 3)]
* Ey[i, j, k])
# Ez component
if nx != 1 or ny != 1:
@@ -170,7 +215,9 @@ cpdef void {{ item.name_b }}(
for k in range(0, nz):
material = ID[2, i, j, k]
for pole in range(maxpoles):
Tz[pole, i, j, k] = Tz[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ez[i, j, k]
Tz[pole, i, j, k] = (Tz[pole, i, j, k]
- updatecoeffsdispersive[material, 2 + (pole * 3)]
* Ez[i, j, k])
{% endfor %}
@@ -199,13 +246,15 @@ cpdef void {{ item.name_a_1 }}(
{{ item.field_type }}[:, :, ::1] Hy,
{{ item.field_type }}[:, :, ::1] Hz
):
"""This function updates the electric field components when dispersive materials (with 1 pole) are present.
"""This function updates the electric field components when dispersive
materials (with 1 pole) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
maxpoles (int): Maximum number of poles
updatecoeffs, T, ID, E, H (memoryviews): Access to update coeffients, temporary, ID and field component arrays
nx, ny, nz: int for grid size in cells.
nthreads: int for number of threads to use.
maxpoles: int for maximum number of poles.
updatecoeffs, T, ID, E, H: memoryviews to access to update coeffients,
temporary, ID and field component arrays.
"""
cdef Py_ssize_t i, j, k
@@ -219,12 +268,21 @@ cpdef void {{ item.name_a_1 }}(
for k in range(1, nz):
material = ID[0, i, j, k]
{% if 'complex' in item.dispersive_type %}
phi = {{ item.real_part }}(updatecoeffsdispersive[material, 0]) * {{ item.real_part }}(Tx[0, i, j, k])
{% if item.iswin %}
phi = (updatecoeffsdispersive[material, 0].real
* Tx[0, i, j, k].real)
{% else %}
phi = ({{ item.real_part }}(updatecoeffsdispersive[material, 0])
* {{ item.real_part }}(Tx[0, i, j, k]))
{% endif %}
{% else %}
phi = updatecoeffsdispersive[material, 0] * Tx[0, i, j, k]
{% endif %}
Tx[0, i, j, k] = updatecoeffsdispersive[material, 1] * Tx[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ex[i, j, k]
Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi
Tx[0, i, j, k] = (updatecoeffsdispersive[material, 1] * Tx[0, i, j, k]
+ updatecoeffsdispersive[material, 2] * Ex[i, j, k])
Ex[i, j, k] = (updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2]
* (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3]
* (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi)
# Ey component
if nx != 1 or nz != 1:
@@ -233,12 +291,21 @@ cpdef void {{ item.name_a_1 }}(
for k in range(1, nz):
material = ID[1, i, j, k]
{% if 'complex' in item.dispersive_type %}
phi = {{ item.real_part }}(updatecoeffsdispersive[material, 0]) * {{ item.real_part }}(Ty[0, i, j, k])
{% if item.iswin %}
phi = (updatecoeffsdispersive[material, 0].real
* Ty[0, i, j, k].real)
{% else %}
phi = ({{ item.real_part }}(updatecoeffsdispersive[material, 0])
* {{ item.real_part }}(Ty[0, i, j, k]))
{% endif %}
{% else %}
phi = updatecoeffsdispersive[material, 0] * Ty[0, i, j, k]
{% endif %}
Ty[0, i, j, k] = updatecoeffsdispersive[material, 1] * Ty[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ey[i, j, k]
Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
Ty[0, i, j, k] = (updatecoeffsdispersive[material, 1] * Ty[0, i, j, k]
+ updatecoeffsdispersive[material, 2] * Ey[i, j, k])
Ey[i, j, k] = (updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3]
* (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1]
* (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi)
# Ez component
if nx != 1 or ny != 1:
@@ -247,12 +314,21 @@ cpdef void {{ item.name_a_1 }}(
for k in range(0, nz):
material = ID[2, i, j, k]
{% if 'complex' in item.dispersive_type %}
phi = {{ item.real_part }}(updatecoeffsdispersive[material, 0]) * {{ item.real_part }}(Tz[0, i, j, k])
{% if item.iswin %}
phi = (updatecoeffsdispersive[material, 0].real
* Tz[0, i, j, k].real)
{% else %}
phi = ({{ item.real_part }}(updatecoeffsdispersive[material, 0])
* {{ item.real_part }}(Tz[0, i, j, k]))
{% endif %}
{% else %}
phi = updatecoeffsdispersive[material, 0] * Tz[0, i, j, k]
{% endif %}
Tz[0, i, j, k] = updatecoeffsdispersive[material, 1] * Tz[0, i, j, k] + updatecoeffsdispersive[material, 2] * Ez[i, j, k]
Ez[i, j, k] = updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi
Tz[0, i, j, k] = (updatecoeffsdispersive[material, 1] * Tz[0, i, j, k]
+ updatecoeffsdispersive[material, 2] * Ez[i, j, k])
Ez[i, j, k] = (updatecoeffsE[material, 0] * Ez[i, j, k] + updatecoeffsE[material, 1]
* (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[material, 2]
* (Hx[i, j, k] - Hx[i, j - 1, k]) - updatecoeffsE[material, 4] * phi)
{% endfor %}
@@ -277,13 +353,15 @@ cpdef void {{ item.name_b_1 }}(
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez
):
"""This function updates a temporary dispersive material array when disperisive materials (with 1 pole) are present.
"""This function updates a temporary dispersive material array when
disperisive materials (with 1 pole) are present.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
maxpoles (int): Maximum number of poles
updatecoeffs, T, ID, E (memoryviews): Access to update coeffients, temporary, ID and field component arrays
nx, ny, nz: int for grid size in cells.
nthreads: int for number of threads to use.
maxpoles: int for maximum number of poles.
updatecoeffs, T, ID, E, H: memoryviews to access to update coeffients,
temporary, ID and field component arrays.
"""
cdef Py_ssize_t i, j, k

41
setup.py
查看文件

@@ -46,7 +46,9 @@ def build_dispersive_material_templates():
functions.
"""
env = Environment(loader = FileSystemLoader('gprMax/templates'), )
iswin = True if sys.platform == 'win32' else False
env = Environment(loader = FileSystemLoader(os.path.join('gprMax', 'templates')), )
template = env.get_template('fields_updates_dispersive_template')
@@ -61,7 +63,8 @@ def build_dispersive_material_templates():
'name_a_1': 'update_electric_dispersive_1pole_A_double_real',
'name_b_1': 'update_electric_dispersive_1pole_B_double_real',
'field_type': 'double',
'dispersive_type': 'double'
'dispersive_type': 'double',
'iswin': iswin
},
# templates for Float precision and dispersive materials with
# real susceptibility functions
@@ -71,7 +74,8 @@ def build_dispersive_material_templates():
'name_a_1': 'update_electric_dispersive_1pole_A_float_real',
'name_b_1': 'update_electric_dispersive_1pole_B_float_real',
'field_type': 'float',
'dispersive_type': 'float'
'dispersive_type': 'float',
'iswin': iswin
},
# templates for Double precision and dispersive materials with
# complex susceptibility functions
@@ -83,7 +87,8 @@ def build_dispersive_material_templates():
'field_type': 'double',
'dispersive_type': 'double complex',
# c function to take real part of complex double type
'real_part': 'creal'
'real_part': 'creal',
'iswin': iswin
},
# templates for Float precision and dispersive materials with
# complex susceptibility functions
@@ -95,15 +100,17 @@ def build_dispersive_material_templates():
'field_type': 'float',
'dispersive_type': 'float complex',
# c function to take real part of complex double type
'real_part': 'crealf'
'real_part': 'crealf',
'iswin': iswin
}]
)
with open('gprMax/cython/fields_updates_dispersive.pyx', 'w') as f:
with open(os.path.join('gprMax', 'cython', 'fields_updates_dispersive.pyx'), 'w') as f:
f.write(r)
# Generate Cython file for dispersive materials update functions
if not os.path.isfile('gprMax/cython/fields_updates_dispersive.pyx'):
cython_disp_file = os.path.join('gprMax', 'cython', 'fields_updates_dispersive.pyx')
if not os.path.isfile(cython_disp_file):
build_dispersive_material_templates()
# Process 'build' command line argument
@@ -143,6 +150,7 @@ if 'cleanall' in sys.argv:
print(f'Removed: {os.path.abspath(libfile)}')
except OSError:
print(f'Could not remove: {os.path.abspath(libfile)}')
# Remove build, dist, egg and __pycache__ directories
shutil.rmtree(Path.cwd().joinpath('build'), ignore_errors=True)
shutil.rmtree(Path.cwd().joinpath('dist'), ignore_errors=True)
@@ -150,6 +158,11 @@ if 'cleanall' in sys.argv:
for p in Path.cwd().rglob('__pycache__'):
shutil.rmtree(p, ignore_errors=True)
print(f'Removed: {p}')
# Remove 'gprMax/cython/fields_updates_dispersive.pyx' if its there
if os.path.isfile(cython_disp_file):
os.remove(cython_disp_file)
print(f'Removed: {cython_disp_file}')
# Now do a normal clean
sys.argv[1] = 'clean' # this is what distutils understands
@@ -210,23 +223,13 @@ else:
extensions = []
for file in cythonfiles:
tmp = os.path.splitext(file)
if tmp[0] == 'gprMax/cython/fields_updates_dispersive' and sys.platform == 'win32':
extension = Extension(tmp[0].replace(os.sep, '.'),
extension = Extension(tmp[0].replace(os.sep, '.'),
[tmp[0] + tmp[1]],
language='c',
include_dirs=[np.get_include()],
extra_compile_args=compile_args,
extra_link_args=linker_args,
libraries=libraries,
define_macros=[("CYTHON_CCOMPLEX", 0)])
else:
extension = Extension(tmp[0].replace(os.sep, '.'),
[tmp[0] + tmp[1]],
language='c',
include_dirs=[np.get_include()],
extra_compile_args=compile_args,
extra_link_args=linker_args,
libraries=libraries)
libraries=libraries)
extensions.append(extension)
# Cythonize - build .c files