gprMax S11 simulations using GPUs

[Murat]

Hi, 

First of all, thank you all for developing gprMax and making it open source. This is great!

I wonder if there is any way to compute antenna s-parameters without using a transmission line source. (for example, via voltage source or others).

Because, gprMax does not support GPU computation when a transmission line source is included in the model, and running complicated antenna structure simulations on CPU takes too long, especially with a smaller spatial discretisation (e.g., dx=0.1 mm, dy=1mm, dz=0.1 mm). Therefore, it would be great if it is possible to calculate antenna parameters using GPUs.

Thanks

[Antonis Giannopoulos]

Hi Murat,

Thanks for your good words.

You can use the GPU and you do not have to use a transmission line source. You need to use a voltage source with some resistance (use the same value as you will have used for the impedance of an ideal transmission line if you want or a 50 Ohm load). You will have to use MATLAB or Python or any scripting language or tool to post process the results but it is easy. There are two ways depending on how well you understand FDTD. One easier than the other.

Common part: You need to use an #rx: output point as otherwise you will not be recording the response you need. The transmission line gives you that as a default due to its formulation. This point should have the same coordinates as the source point. You should use the #voltage_source: instead of the #transmission_line: 

Example modification of the antenna_wire_dipole_fs.in

#waveform: gaussian 1 1e9 mypulse
#voltage_source: z 0.025 0.025 0.100 73 mypulse
#rx: 0.025 0.025 0.100

1st Approach 

Use as your #waveform: the same gaussian excitation as used in the example (see above). In this case once you import your correct electric field component into your processing script you will need to multiply it by -1 x DL where DL is the correct spatial resolution. This will convert it to Voltage from Electric field at this point. This is the total response but for S11 you need the incident response and the reflected. You need to generate the correct incident waveform (it is very easy and there are description and tools available) sampled at the same DT points as in the FDTD solution. Once you have this one then create the reflected response by subtracting the total voltage response you have got from gprMax the incident one you have calculated (the gaussian pulse). You then FFT the incident and the reflected responses and proceed to calculate S11 as |V_ref|/|V_inc| your solution should be valid for the frequencies your model is valid but you need to restrict your plotting to them. FDTD has a range of validity that depends on the spatial step in relation to the wavelength. The solution is not correct for all frequencies up to Nyquist! 

2nd Approach (not documented !!!)

DISCLAIMER: if you do not understand what you are doing this might confuse you but it is easier to use ...

Use as your #waveform: the description impulse (it is not in the manual) You must use the other parameters in the definition (amplitude is 1 and frequency any number as it is ignored). Once you import your correct electric field component into your processing script you will need to multiply it by -1 x DL where DL is the correct spatial resolution. This will convert it to Voltage from Electric field at this point. This is the total response but for S11 you need the incident response and the reflected. However, the incident is only 1 followed by zeros as it is an impulse. Further, due to the way gprMax records the output this application of the impulse is not recorded and computationally the obtained waveform is directly the reflected response due to such an impulse excitation. In essence, it is the modelled impulse response of the antenna. This is the FDTD impulse response and not the theoretical one! Looking at this in the time domain it makes no much sense. If you convolve it with a proper excitation (e.g a gaussian) it will but you do not need to do this to get S11. You just need to FFT this response and your S11 is this response! As S11 is |V_ref|/|V_inc| and on this case |V_inc| = 1. You need to be careful as this is not valid for all frequencies and it will look strange and wrong if you plot them all!. It is valid up to the frequency that your specific FDTD model resolution allows! Be warned!

Example modification of the antenna_wire_dipole_fs.in

#waveform: impulse 1 1e9 mypulse
#voltage_source: z 0.025 0.025 0.100 73 mypulse 
#rx: 0.025 0.025 0.100

Hope this works for you and the speed up using the GPU will make this fast enough for your work!

Best and #besafe

Antonis

[Murat]

Hi Antonis,

Thank you very much for explaining these two methods in details.

 

These solutions sound quite logical and easy to implement. I will simulate the antennas and compute their S11 parameters using these methods and also compare the results to the ones obtained via transmission line source.

Best wishes,

Murat

[Murat]

Hi Antonis,

I have tried the methods you suggested but there is a difference between s11 parameters obtained using transmission line and voltage as shown in the figure below.  I used a bowtie antenna model for comparison.  I tried many times but could not figure out why they are different. I have the codes below for transmission line and voltage source models and processing them. Could you please let me know if there is anything I have missed during simulations or calculations.

Thank you very much.

Figure: S11 [dB] vs Frequency [Hz] - Blue line indicates the S11 calculated using a voltage source (impulse waveform) and the orange line indicates the S-11 calculated using a transmission line (impulse waveform).    

Codes:

#Waveform and sources

print('#waveform: impulse 1 1.5e9 mypulse')

#print('#transmission_line: x {:g} {:g} {:g} 50 mypulse'.format(tx_pos[0], tx_pos[1], tx_pos[2]))  # for comparison

print('#voltage_source: x {:g} {:g} {:g} 50 mypulse'.format(tx_pos[0], tx_pos[1], tx_pos[2]))

rx(tx_pos[0], tx_pos[1], tx_pos[2])

S11 calculations with Transmission Line Source:

In the transmission line model, we can see the incident voltage from the output data as indicated Vinc below. However, When I checked it, it is not an impulse, I guess this is caused by the transmission line, therefore, I subtracted it from the total wave received.

rf = h5py.File('./Bowtie.out', 'r')

# read the spatial resolutions of the model and simulations 

dxdydz = rf.attrs['dx_dy_dz']

# read the time resolution

dt = rf.attrs['dt']

# read number of iterations

iterations = rf.attrs['Iterations']

# simulation time array

time = np.linspace(0, (iterations - 1) * dt, num=iterations)

# simulation frequency resolution

df = 1 / np.amax(time)

#rx path in the results file

rxpath = '/rxs/rx1/'

#get E-Fields

Ex = rf[rxpath + 'Ex'][:]

Ey = rf[rxpath + 'Ey'][:]

Ez = rf[rxpath + 'Ez'][:]

#calculate Vx Vy Vz

Vx = -1*Ex*dxdydz[0]

Vy = -1*Ey*dxdydz[1]

Vz = -1*Ez*dxdydz[2]

# use the antenna parameter calculation script of gprMax

antennaparams = calculate_antenna_params('./Bowtie.out') 

 

Vtotal = antennaparams['Vtotal']   # Output using TL model

Vinc = antennaparams['Vinc']   # Incident wave

Vref = Vx - Vinc     #    Vref =  Vtotal - Vinc   also gives the same result via the antenna parameter calculation script of gprMax

freqs = np.fft.fftfreq(Vinc.size, d=dt)

with np.errstate(divide='ignore'):

        s11 = np.abs(np.fft.fft(Vref) / np.fft.fft(Vinc))

s11 = 20 * np.log10(s11)

plt.plot(freqs[0:200], s11[0:200])

S11 calculations with Voltage Source:  

rf = h5py.File('./Bowtie.out', 'r')

# read the spatial resolutions of the model and simulations 

dxdydz = rf.attrs['dx_dy_dz']

# read the time resolution

dt = rf.attrs['dt']

# read number of iterations

iterations = rf.attrs['Iterations']

# simulation time array

time = np.linspace(0, (iterations - 1) * dt, num=iterations)

# simulation frequency resolution

df = 1 / np.amax(time)

#rx path in the results file

rxpath = '/rxs/rx1/'

#get E-fields

Ex = rf[rxpath + 'Ex'][:]

Ey = rf[rxpath + 'Ey'][:]

Ez = rf[rxpath + 'Ez'][:]

#calculate Vx, Vy, Vz

Vx = -1*Ex*dxdydz[0]

Vy = -1*Ey*dxdydz[1]

Vz = -1*Ez*dxdydz[2]

#impulse is not recorded in the Vx, so we can use Vx directly to calculate the S11

fftVx = np.fft.fft(Vx)

s11x = 20 * np.log10(fftVx)

freqs = np.fft.fftfreq(Vx.size, d=dt)

plt.plot(freqs[0:200], s11x[0:200])

Note:

1. In transmission line model, Vx (obtained via Rx with Ex-field) and Vtotal (obtained through the transmission line)  are exactly the same so I believe I do not have any problem with obtaining the received/total wave voltage.

2. Although the transmission line is fed by an impulse, Vinc is not impulse as shown in the figure below. So, I guess incident wave using voltage source might also be different.

Best regards,

Murat

[Antonis Giannopoulos]

Hi Murat,

Interesting results. We have not used the transmission line with an impulse but I don't think in this case should have been a problem. Maybe running the transmission line with the gaussian and comparing the results you get with the impulse it might be a good place to see if the two methods differ in something that have been missed. The impulse response in the time domain will look strange as the FDTD model does not properly support impulsive propagation. However, when you convolve this response with a proper forcing function that has an acceptable frequency content, for the given model, the result should be fine as in essence you filter out the erroneous parts. In a similar way the FFT part of this impulse response should be OK for the "valid" frequencies and obviously not for frequencies that are not properly resolved.

The other thing to check is to use the voltage source with a gaussian excitation but you will need to properly recreate the incident field in Python to do the calculations. This should be easy. If this compares well with the transmission line then we need to look a bit more closely on how the impulse does not do exactly what is needed. 

Best

Antonis

[Murat]

Hi Antonis 

Thank you very much for your suggestions and help.

I simulated only a transmission line with a gaussian input wave and one Rx at the same location without having any structure to see if Vinc (calculated from transmission line) and Vx   (=-1*Ex*dx) are the same. When I compared them, I realized that Vinc = Vx/2.  

Subsequently, I simulated only a voltage source and an Rx at the same location without having any structure and recorded the waveform which is denoted by VSinc (=-1*Ex*dx) below in the codes.  During antenna simulations using voltage source and Rx, I used VSinc/2 as the incident wave to calculate S11, it seemed that the results were correct as S11 parameters obtained through this method and using a transmission line were quite similar as shown in the figure below. 

Many thanks,

Murat

VSinc = np.loadtxt('./VS_Vinput.txt')     # Gaussian waveform obtained from the simulation consisting of only a Voltage Source and an Rx (Vmax =1)

TLS11 = np.loadtxt('./TLS11.txt')           #s11 results of the antenna obtained through the simulation via a transmission line 

# Calculation of S11 of the antenna simulated via a voltage source and Rx

freqs = np.fft.fftfreq(Vx.size, d=dt)

Vref = Vx-VSinc/2

S11 = np.fft.fft(Vref)/np.fft.fft(VSinc/2)

S11dB = 20 * np.log10(S11)

plt.plot(freqs[0:200],S11dB[0:200])

plt.plot(freqs[0:200],TLS11[0:200])

 

[Antonis Giannopoulos]

Hi Murat,

That makes some sense as looking down the line from a Thevenin equivalent point of view you can replace it with a 2V and an admittance across. It is good to try and more detailed comparisons to find some equivalent sourcing models. There is a paper about this that is well known and useful 

10.1109/TAP.2003.815414

The idea of using a 1D transmission line model for such excitation is not as ideal as it is when used as a model for a coax cable feed on a monopole but it is kind of OK. Results should be equivalent but often small phase variations can exist depending on different excitations.

Best

Antonis

[Murat]

Hi Antonis,

Thanks for the explanation and the paper which seems quite useful for understanding FDTD feed models.  I will look into this more by considering various antenna types.

Best,

Murat