diff --git a/gprMax/templates/fields_updates_dispersive_template b/gprMax/templates/fields_updates_dispersive_template index ec9e5c58..8c49d69c 100644 --- a/gprMax/templates/fields_updates_dispersive_template +++ b/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 diff --git a/setup.py b/setup.py index 9eded61b..f905b687 100644 --- a/setup.py +++ b/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