Question on #excitation_file (excitation file)

[Fabio da Silva]

Good morning,

I have another question, this time about excitation files. I have a waveform I want to generate comprising a sine wave enveloped by a Gaussian pulse. The reasons for this choice rest on retaining as much spectral purity of the sine wave as possible.

Here is where the problem lies. When I try to execute the code, I get:

For the pulse containing sine and envelope:

FloatingPointError: divide by zero encountered in log10

Now, if I simply take the one component off the pulse (sine or envelope) and send the other function only, it works with the following warning issued:

WARNING: Numerical dispersion analysis not carried out as either no waveform detected or waveform does not fit within specified time window and is therefore being truncated.

Input and output files are attached for your convenience.

Disclaimer: because I am testing the boundaries of my hardware, there is a large number of cells in the simulation, so you may need to shrink it for speed.

Thanks for any insight on this,

Fabio.

PS: I noticed that I get the FloatingPointError for all three wave forms when I try to make the excitation file multi column.

[Antonis Giannopoulos]

Hi, Fabio

A user pulse needs to be sampled at the correct time instance (actually this is at half step intervals (i.e. 0.5*Dt, 1.5*Dt, ...) if this is a Hertzian dipole. Further, it should have a length equal to the time window (2000 iterations) or otherwise is going to be truncated with zeros. In your case, this is not going to give much tone purity as your pulse is not going to zero. I think you need to sample for more so no visible ripples are observed when the pulse stops. This will improve the purity. I am not sure why the frequency dispersion determination function gives an error but I suppose Craig should be able to look at this.

Best

Antonis

[Fabio da Silva]

Thanks, Antonis.

You are correct in your assessment of tone purity and how the pulse I used is not pure according to my goals. In the real simulation, your observation will surely be implemented. I was not aware that the time evolution evaluations occur at half time steps. Still, since their separation is the same as the full time steps, all I need to do is make sure I add a half time step offset to the time analysis and I will be aligned. Incidentally, I use a time interval specified by the number of time steps and not actual time values to make sure I use the internal definitions of gprMax.

Back to the numerical dispersion problem, I can surely live with it as long as I understand the logic behind its exception call. I tried to glean it from different functions I used, but could not figure it out. In fact the pulse example I was the best way I found to explain the problem. I am glad to add any further information you need.

Cordially,

Fabio.

[Fabio da Silva]

The exception happens due to a division by zero inside the log10 function. I did a grep into the gprMax folder for "log10" and got the following (the results with a division inside the log10 call are highlighted in yellow:

C:\Users\fcss\Documents\gprMax>grep -R log10 *

build/lib.win-amd64-3.5/gprMax/grid.py:                    power = 10 * np.log10(np.abs(np.fft.fft(waveformvalues))**2)

build/lib.win-amd64-3.5/tests/test_models.py:            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])  # Ignore any zero division in log10

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    Vincp = 20 * np.log10(np.abs((np.fft.fft(Vinc) * delaycorrection)))

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    Iincp = 20 * np.log10(np.abs(np.fft.fft(Iinc)))

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    Vrefp = 20 * np.log10(np.abs((np.fft.fft(Vref) * delaycorrection)))

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    Irefp = 20 * np.log10(np.abs(np.fft.fft(Iref)))

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    Vtotalp = 20 * np.log10(np.abs((np.fft.fft(Vtotal) * delaycorrection)))

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    Itotalp = 20 * np.log10(np.abs(np.fft.fft(Itotal)))

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:    s11 = 20 * np.log10(s11)

build/lib.win-amd64-3.5/tools/plot_antenna_params.py:        s21 = 20 * np.log10(s21)

build/lib.win-amd64-3.5/tools/plot_Ascan.py:                power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2)

build/lib.win-amd64-3.5/tools/plot_source_wave.py:        start = np.where((10 * np.log10(waveform / np.amax(waveform))) > powerdrop)[0][0]

build/lib.win-amd64-3.5/tools/plot_source_wave.py:        stop = np.where((10 * np.log10(waveform[start:] / np.amax(waveform))) < powerdrop)[0][0] + start

build/lib.win-amd64-3.5/tools/plot_source_wave.py:        power = 10 * np.log10(np.abs(np.fft.fft(waveform))**2)

build/lib.win-amd64-3.5/user_libs/antenna_patterns/plot_fields.py:    ax.plot(theta, 10 * np.log10(pattplot), label='{:.2f}m'.format(radii[patt]), marker='.', ms=6, lw=1.5)

build/lib.win-amd64-3.5/user_libs/antenna_patterns/plot_fields.py:# ax.plot(theta, 10 * np.log10(hertzplot1), label='Inf. dipole, 0.1m', color='black', ls='-.', lw=3)

build/lib.win-amd64-3.5/user_libs/antenna_patterns/plot_fields.py:# ax.plot(theta, 10 * np.log10(hertzplot2), label='Inf. dipole, 0.58m', color='black', ls='--', lw=3)

build/lib.win-amd64-3.5/user_libs/optimisation_taguchi/fitness_functions.py:            tmp = 20 * np.log10(np.abs(modelresp - refresp) / np.amax(np.abs(refresp)))

gprMax/grid.py:                    power = 10 * np.log10(np.abs(np.fft.fft(waveformvalues))**2)

Binary file gprMax/__pycache__/grid.cpython-35.pyc matches

tests/test_models.py:            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])  # Ignore any zero division in log10

tools/plot_antenna_params.py:    Vincp = 20 * np.log10(np.abs((np.fft.fft(Vinc) * delaycorrection)))

tools/plot_antenna_params.py:    Iincp = 20 * np.log10(np.abs(np.fft.fft(Iinc)))

tools/plot_antenna_params.py:    Vrefp = 20 * np.log10(np.abs((np.fft.fft(Vref) * delaycorrection)))

tools/plot_antenna_params.py:    Irefp = 20 * np.log10(np.abs(np.fft.fft(Iref)))

tools/plot_antenna_params.py:    Vtotalp = 20 * np.log10(np.abs((np.fft.fft(Vtotal) * delaycorrection)))

tools/plot_antenna_params.py:    Itotalp = 20 * np.log10(np.abs(np.fft.fft(Itotal)))

tools/plot_antenna_params.py:    s11 = 20 * np.log10(s11)

tools/plot_antenna_params.py:        s21 = 20 * np.log10(s21)

tools/plot_Ascan.py:                power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2)

tools/plot_source_wave.py:        start = np.where((10 * np.log10(waveform / np.amax(waveform))) > powerdrop)[0][0]

tools/plot_source_wave.py:        stop = np.where((10 * np.log10(waveform[start:] / np.amax(waveform))) < powerdrop)[0][0] + start

tools/plot_source_wave.py:        power = 10 * np.log10(np.abs(np.fft.fft(waveform))**2)

Binary file tools/__pycache__/plot_Ascan.cpython-35.pyc matches

user_libs/antenna_patterns/plot_fields.py:    ax.plot(theta, 10 * np.log10(pattplot), label='{:.2f}m'.format(radii[patt]), marker='.', ms=6, lw=1.5)

user_libs/antenna_patterns/plot_fields.py:# ax.plot(theta, 10 * np.log10(hertzplot1), label='Inf. dipole, 0.1m', color='black', ls='-.', lw=3)

user_libs/antenna_patterns/plot_fields.py:# ax.plot(theta, 10 * np.log10(hertzplot2), label='Inf. dipole, 0.58m', color='black', ls='--', lw=3)

user_libs/optimisation_taguchi/fitness_functions.py:            tmp = 20 * np.log10(np.abs(modelresp - refresp) / np.amax(np.abs(refresp)))

user_models/sd.exe:extern __declspec(__host__) __declspec(__device__) __declspec(__device_builtin__) double         __cdecl log10(double x) ;

user_models/sd.exe:extern  __declspec(__host__) __declspec(__device__) __declspec(__device_builtin__) float __cdecl log10f(float) ;

user_models/sd.exe:      double __cdecl log10(  double _X);

user_models/sd.exe:           float  __cdecl log10f(  float _X);

user_models/sd.exe:      __inline long double __cdecl log10l(  long double _X)

user_models/sd.exe:        return log10((double)_X);

user_models/sd.exe:  inline float log10(  float _Xx) noexcept

user_models/sd.exe:     return (:: log10f(_Xx));

user_models/sd.exe:  inline long double log10(  long double _Xx) noexcept

user_models/sd.exe:     return (:: log10l(_Xx));

user_models/sd.exe:extern "C"    double __cdecl log10(  double); template<class _Ty> inline typename ::std:: enable_if< ::std:: is_integral<_Ty>::value, double>::type log10(_Ty _Left) { return (:: log10((double)_Left)); }

user_models/sd.exe:using :: log10; using :: modf; using :: pow;

user_models/sd.exe:using :: log10f; using :: modff; using :: powf;

user_models/sd.exe:using :: log10l; using :: modfl; using :: powl;

user_models/sd.exe:extern __declspec(__host__) __declspec(__device__) __declspec(__cudart_builtin__) float    __cdecl log10(float) throw();

user_models/sd.exe:__declspec(__device__) __declspec(__cudart_builtin__) __declspec(__device_builtin__) __declspec(__cudart_builtin__) float                  __log10f(float x) ;

user_models/sd.exe:__declspec(__device__) complex<float> log10(const complex<float>&);

user_models/sd.exe:__declspec(__device__) complex<double> log10(const complex<double>&);

I am still studying the code and will report any findings as soon as I get them.

Fabio.

[Craig Warren]

Hi Fabio,

The log10 evaluated on zero error you are experiencing comes from the dispersion analysis function (in the grid.py module if you are interested). The function is there to help guide users as to whether their model will likely suffer from the effects of numerical dispersion. I have put in a fix for the issue to catch the case when there is a zero from the FFT. Anyway if you update gprMax (pull from GitHub) the error should now be fixed.

Kind regards,

Craig

[Fabio da Silva]

Thanks, Craig.

It worked like a charm. I really appreciate your prompt response.

Cordially,

Fabio.