From 477e5c4fe296f4740eb721810a4ae10b1170e60b Mon Sep 17 00:00:00 2001 From: AdittyaPal <42172854+AdittyaPal@users.noreply.github.com> Date: Sun, 13 Aug 2023 19:37:31 +0530 Subject: [PATCH] Updated the (main-like) function --- gprMax/cython/planeWaveModules.pyx | 293 ++++++++++++++++++++++++----- 1 file changed, 249 insertions(+), 44 deletions(-) diff --git a/gprMax/cython/planeWaveModules.pyx b/gprMax/cython/planeWaveModules.pyx index 972ac2a3..25848326 100644 --- a/gprMax/cython/planeWaveModules.pyx +++ b/gprMax/cython/planeWaveModules.pyx @@ -1,7 +1,9 @@ + 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 @cython.cdivision(True) @@ -101,6 +103,21 @@ cdef inline double getTanValue(int[:] integers, double[:] dr): @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)], @@ -119,10 +136,10 @@ cdef int[:] getTheta(int m_x, int m_y, double theta, double Delta_theta, double[ 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 : projections of the propagation vector along the different coordinate axes + 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: @@ -169,7 +186,8 @@ def getIntegerForAngles(double phi, double Delta_phi, double theta, double Delta Returns: -------------------------- - integers, int array : sequence of the three integers [m_x, m_y, m_z] + directions, int array : the integers specifying the direction of propagation of the plane wave along the three coordinate axes + integers, 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 @@ -210,6 +228,25 @@ def getIntegerForAngles(double phi, double Delta_phi, double theta, double Delta @cython.boundscheck(False) cdef double[:, :, :, :, ::1] updateHFields(double[:] coefficients, int n_x, int n_y, int n_z, double[:, :, :, :, ::1] fields, int waveID): + ''' + Method to update the magnetic fields in the TFSF grid uding the electric fields from the previous update. + __________________________ + + Input parameters: + -------------------------- + coefficients, float : stores the coefficients of the fields in the update equation for the magnetic field + n_x, int : stores the number of grid cells along the x axis of the TFSF box + n_y, int : stores the number of grid cells along the y axis of the TFSF box + n_z, int : stores the number of grid cells along the z axis of the TFSF box + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + waveID, int : stores the index number of the field in case of multiple plane waves + __________________________ + + Returns: + -------------------------- + fields, double array : returns the updated magnetic fields for the grid cells over the TFSF box at particular indices + ''' cdef Py_ssize_t i, j, k = 0 cdef double c0 = coefficients[0] @@ -245,6 +282,25 @@ cdef double[:, :, :, :, ::1] updateHFields(double[:] coefficients, int n_x, int @cython.boundscheck(False) cdef double[:, :, :, :, ::1] updateEFields(double[:] coefficients, int n_x, int n_y, int n_z, double[:, :, :, :, ::1] fields, int waveID): + ''' + Method to update the electric fields in the TFSF grid uding the electric fields from the previous update. + __________________________ + + Input parameters: + -------------------------- + coefficients, float : stores the coefficients of the fields in the update equation for the electric field + n_x, int : stores the number of grid cells along the x axis of the TFSF box + n_y, int : stores the number of grid cells along the y axis of the TFSF box + n_z, int : stores the number of grid cells along the z axis of the TFSF box + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + waveID, int : stores the index number of the field in case of multiple plane waves + __________________________ + + Returns: + -------------------------- + fields, double array : returns the updated electric fields for the grid cells over the TFSF box at particular indices + ''' cdef Py_ssize_t i, j, k = 0 cdef double c0 = coefficients[0] @@ -285,6 +341,27 @@ cdef double[:, :, :, :, ::1] updateEFields(double[:] coefficients, int n_x, int @cython.boundscheck(False) cdef double[:, :, :, ::1] implementABC(double[:, :, :, ::1] face_fields, double abccoef, int n_x, int n_y, int n_z, double[:, :, :, :, ::1] fields, int waveID): + ''' + Method to implement the absorbing boundary conditions to the FDTD grid. + __________________________ + + Input parameters: + -------------------------- + face_fields, double array : stores the fields on the six faces of the FDTD Box + abccoef, double : stores the absorbing boundary coefficient so that there are no reflection from the edge of the grid + n_x, int : stores the number of grid cells along the x axis of the TFSF box + n_y, int : stores the number of grid cells along the y axis of the TFSF box + n_z, int : stores the number of grid cells along the z axis of the TFSF box + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + waveID, int : stores the index number of the field in case of multiple plane waves + __________________________ + + Returns: + -------------------------- + face_fields, double array : returns the updated fields on the surfaces of the FDTD Box + + ''' cdef Py_ssize_t i, j, k = 0 cdef double[:, :, ::1] e_x = fields[waveID, 0, :, :, :] @@ -357,6 +434,27 @@ cdef double[:, :, :, ::1] implementABC(double[:, :, :, ::1] face_fields, double @cython.boundscheck(False) cdef double[:, :, :, :, ::1] applyTFSFMagnetic(double[:] coefficients, double[:, :] e1D, int[:] m, double[:, :, :, :, ::1]fields, int[:, :] corners, int waveID): + ''' + 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: + -------------------------- + coefficients, float : stores the coefficients of the fields in the TFSF assignment equation for the magnetic field + e1D, 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 + 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) + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries + waveID, int : stores the index number of the field in case of multiple plane waves + __________________________ + + Returns: + -------------------------- + fields, double array : returns the updated fields atre applying the TF/Sf formulation at the boundary of the TF/SF surface + + ''' cdef Py_ssize_t i, j, k = 0 cdef int x_start = corners[0, 0] @@ -462,6 +560,27 @@ cdef double[:, :, :, :, ::1] applyTFSFMagnetic(double[:] coefficients, double[:, cdef double[:, :, :, :, ::1] applyTFSFElectric(double[:] coefficients, double[:, :] h1D, int[:] m, double[:, :, :, :, ::1] fields, int[:, :] corners, int waveID): + ''' + 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: + -------------------------- + coefficients, float : stores the coefficients of the fields in the TFSF assignment equation for the electric field + h1D, double array : array to store the magnetic fields of the one dimensional representation of the plane wave in a direction along which the wave propagates + 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) + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries + waveID, int : stores the index number of the field in case of multiple plane waves + __________________________ + + Returns: + -------------------------- + fields, double array : returns the updated fields atre applying the TF/Sf formulation at the boundary of the TF/SF surface + + ''' cdef Py_ssize_t i, j, k = 0 cdef int x_start = corners[0, 0] cdef int y_start = corners[0, 1] @@ -561,7 +680,13 @@ cdef double[:, :, :, :, ::1] applyTFSFElectric(double[:] coefficients, double[:, e_x[i, j, k] -= coefficients[2] * h1D[1, index] return fields -''' + +@cython.cdivision(True) +@cython.wraparound(False) +@cython.boundscheck(False) +cdef double[:, ::1] updateMagneticFields(int n, double[:] H_coefficients, double[:, ::1] H_fields, + double[:, ::1] E_fields, int dimensions, int[:] m): + ''' Method to update the magnetic fields for the next time step using Equation 8 of DOI: 10.1109/LAWP.2009.2016851 __________________________ @@ -572,20 +697,17 @@ cdef double[:, :, :, :, ::1] applyTFSFElectric(double[:] coefficients, double[:, 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 + dimensions, int : stores the number of dimensions in which the simulation is run + 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: -------------------------- H_fields, double array : magnetic field array with the axis entry for the current time added -''' -@cython.cdivision(True) -@cython.wraparound(False) -@cython.boundscheck(False) -cdef double[:, ::1] updateMagneticFields(int n, double[:] H_coefficients, double[:, ::1] H_fields, - double[:, ::1] E_fields, int dimensions, int[:] m): - + ''' cdef Py_ssize_t i, j = 0 cdef int dim_mod1, dim_mod2 = 0 for i in prange(dimensions, nogil=True, schedule='static'): #loop to update each component of the magnetic field @@ -613,7 +735,10 @@ cdef double[:, ::1] updateElectricFields(int n, double[:] E_coefficients, double 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 + dimensions, int : stores the number of dimensions in which the simulation is run + 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: @@ -633,29 +758,81 @@ cdef double[:, ::1] updateElectricFields(int n, double[:] E_coefficients, double return E_fields @cython.cdivision(True) -cdef inline double getSource(double time, double ppw): +cdef double getSource(double time, double ppw, 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 - ppw, int : points per wavelength for the wave source + time, float : time at which the magnitude of the source is calculated + ppw, int : 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 array : magnitude of the source for the requested indices at the current time + sourceMagnitude, double : magnitude of the source for the requested indices at the current time ''' - return (1.0-2.0*(M_PI*(time/ppw - 1.0))*(M_PI*(time/ppw - 1.0)) - )*exp(-(M_PI*(time/ppw - 1.0))*(M_PI*(time/ppw - 1.0))) # define a Ricker wave source + # Waveforms + if (strcmp(wavetype, "gaussian") == 0): + return exp(-(M_PI*(time/ppw - 1.0)) * (M_PI*(time/ppw - 1.0))) + + elif (strcmp(wavetype, "gaussiandot") == 0 or strcmp(wavetype, "gaussianprime") == 0): + return 4 * M_PI/ppw * (M_PI*(time/ppw - 1.0) + ) * exp(-(M_PI*(time/ppw - 1.0)) * (M_PI*(time/ppw - 1.0))) + + elif (strcmp(wavetype, "gaussiandotnorm") == 0): + return 2 * (M_PI*(time/ppw - 1.0) + ) * exp(-(M_PI*(time/ppw - 1.0)) * (M_PI*(time/ppw - 1.0))) + + elif (strcmp(wavetype, "gaussiandotdot") == 0 or strcmp(wavetype, "gaussiandoubleprime") == 0): + return 2 * M_PI*M_PI/(ppw*ppw) * (1.0-2.0*(M_PI*(time/ppw - 1.0))*(M_PI*(time/ppw - 1.0)) + ) * exp(-(M_PI*(time/ppw - 1.0)) * (M_PI*(time/ppw - 1.0))) + + elif (strcmp(wavetype, "gaussiandotdotnorm") == 0 or strcmp(wavetype, "ricker") == 0): + return (1.0 - 2.0*(M_PI*(time/ppw - 1.0)) * (M_PI*(time/ppw - 1.0)) + ) * exp(-(M_PI*(time/ppw - 1.0)) * (M_PI*(time/ppw - 1.0))) # define a Ricker wave source + + elif (strcmp(wavetype, "sine") == 0): + return sin(2 * M_PI * time/ppw) + + elif (strcmp(wavetype, "contsine") == 0): + return min(0.25 * time/ppw, 1) * sin(2 * M_PI * time/ppw) + + elif (strcmp(wavetype, "impulse") == 0): + if (time < dt): # time < dt condition required to do impulsive magnetic dipole + return 1 + else: + return 0 @cython.wraparound(False) @cython.boundscheck(False) cdef double[:, :, :, :, ::1] superImposeFields(planeWaves, int n_x, int n_y, int n_z, int dimensions, double[:, :, :, :, ::1] fields, int iterate): + ''' + Method to superimpose multiple plane waves if the user requires multiple plane wave incidences. + __________________________ + + Input parameters: + -------------------------- + planeWaves, list : list of class DiscretePlaneWave() objects, each representing a particular plane wave + n_x, int : stores the number of grid cells along the x axis of the TFSF box + n_y, int : stores the number of grid cells along the y axis of the TFSF box + n_z, int : stores the number of grid cells along the z axis of the TFSF box + dimensions, int : stores the number of dimensions in which the simulation is run + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + iterate, int : stores the number of incident plane waves in case of multiple plane waves + __________________________ + + Returns: + -------------------------- + fields, double array : returns the superimposed fields in the TFSF box as the fields in the last index. + + ''' cdef Py_ssize_t i, j, k, l, n = 0 cdef int x_start, y_start, z_start for n in range(iterate): @@ -675,45 +852,71 @@ cdef double[:, :, :, :, ::1] superImposeFields(planeWaves, int n_x, int n_y, int fields[iterate, l, i, j, k] += fields[n, l, abs(x_start-i), abs(y_start-j), abs(z_start-k)] return fields -''' -@cython.wraparound(False) -@cython.boundscheck(False) -cdef void takeSnapshot3D(double[:, :, ::1] field, char* filename): - cdef FILE *fptr = fopen(filename, "wb") - cdef Py_ssize_t n_x = field.shape[0] - cdef double[:, ::1] slice_view - cdef Py_ssize_t i = 0 - if fptr is NULL: - raise ValueError("Failed to open the file.") - else: - for i in range(n_x): - slice_view = field[i, ...] - fwrite(&slice_view[0, 0], sizeof(double), slice_view.size, fptr) - fclose(fptr) -''' + + @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) + def getGridFields(planeWave, double[:] C, double[:] D, int snapshot, int n_x, int n_y, int n_z, double[:, :, :, :, ::1] fields, int[:, ::1] corners, int time_duration, double[:, :, :, ::1] face_fields, - double abccoef, double dt, int noOfWaves, double c, double ppw, int dimensions): - + double abccoef, double dt, int noOfWaves, double c, double ppw, + int dimensions, filepath, snapFieldNumber, snapFieldSteps): + ''' + Method to implement the time stepping loop of the FDTD simulation + __________________________ + + Input parameters: + -------------------------- + planeWave, list : list of class DiscretePlaneWave() objects, each representing a particular plane wave + 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 + snapshot, int : stores the interval at which a snapshot of the fields is stored + n_x, int : stores the number of grid cells along the x axis of the TFSF box + n_y, int : stores the number of grid cells along the y axis of the TFSF box + n_z, int : stores the number of grid cells along the z axis of the TFSF box + fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order + E_x, E_y, E_z, H_x, H_y, H_z + corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries + time_dimension, int : stores the number of time steps over which the FDTD simulation is run + face_fields, double array : stores the fields on the six faces of the FDTD Box + abccoef, double : stores the absorbing boundary coefficient so that there are no reflection from the edge of the grid + noOfWaves, int : stores the number of plane waves to be incident on the TFSF Box + c, double : speed of light + ppw, double : stores the number of points per wavelength for the requested source + dimensions, int : stores the number of dimensions in which the simulation is run + filepath, string : specifies the file location where the fields are to be written + snapFieldNumber, list: list to store which fields are to be written in the file + snapFieldSteps, list : list to store any specific time steps at which snapshots are to be taken. + + ''' cdef double real_time = 0.0 # initialize the real time to zero # Py_ssize_t is the proper C type for Python array indices - cdef Py_ssize_t i, j, r = 0 #E and H are the one-dimensional waves + cdef Py_ssize_t i, j, k, l, m = 0 #E and H are the one-dimensional waves for i in range(time_duration): for j in range(noOfWaves): for k in range(dimensions): for l in range(planeWave[j].m[dimensions]): #loop to assign the source values of magnetic field to the first few gridpoints planeWave[j].H_fields[k, l] = planeWave[j].projections[k]*getSource(real_time - (l+( planeWave[j].m[(k+1)%dimensions]+planeWave[j].m[(k+2)%dimensions] - )*0.5)*planeWave[j].ds/c, ppw) + )*0.5)*planeWave[j].ds/c, ppw, + planeWave[j].waveformID, planeWave[j].dt) planeWave[j].H_fields = updateMagneticFields(planeWave[j].length, D, planeWave[j].H_fields, planeWave[j].E_fields, dimensions, planeWave[j].m) #update the magnetic fields in the rest of the DPW grid @@ -730,8 +933,10 @@ def getGridFields(planeWave, double[:] C, double[:] D, int snapshot, real_time += dt #take a half time step because the magnetic field and electric field grids are staggered in time - if (i % snapshot == 0): + if (i % snapshot == 0 or i in snapFieldSteps): fields = superImposeFields(planeWave, n_x, n_y, n_z, dimensions, np.asarray(fields), noOfWaves) - takeSnapshot3D(fields[noOfWaves, 0, :, :, :], f'./snapshots/electric_2{i}.dat'.encode('UTF-8')) - fields[noOfWaves, :, :, :, :] = 0 + for m in snapFieldNumber: + takeSnapshot3D(fields[noOfWaves, m, :, :, :], (f'{filepath}_{m}{i}.dat').encode('UTF-8')) + + fields[noOfWaves, :, :, :, :] = 0 #erase the cells bearing the last index so that it does not interfere with the superimposiion fields at some other time instance