How to use PEC/PMC boundary mirroring (Python)?

[VolkerMuehlhaus]

Hello forum,

I would like to try NF2FF boundary mirroring, which is available according to this documentation:
https://docs.openems.de/python/openEMS/nf2ff.html

But I am struggling with the actual use in my Python model. Does anyone have any code example how to use PEC or PMC mirroring?

Best regards
Volker

[thliebig]

Looking at the code I think you do not need to do anything special. If you use e.g. a PMC boundary somewhere the nf2ff will detect it and set mirroring accordingly...
But I'm not sure if I truly ever tried it or what happens if you use PEC e.g. as gnd plane.

[VolkerMuehlhaus]

Hi Thorsten,

many thanks for looking into this!

I tried to implement mirroring using a PEC boundary in the center plane of a yagi antenna. Which also means that the lumped port is between the drawn radiator and the PEC boundary.

I had tested that symmetry model idea using PEC in Empire EM simulation (worked fine), but it fails with my openEMS model below. The voltage blows up to infinity in simulation, so results are nonsense. I think there is something missing, or I made a mistake.

The "full" antenna model (without using symmetry) had worked fine, the issue is only with the "cut-in-half" symmetry version.

Best regards
Volker

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

### Import Libraries
import os, tempfile, math, shutil
from pylab import *

from CSXCAD import CSXCAD

from openEMS import openEMS
from openEMS.physical_constants import *

#
# FUNCTION TO GIVE RANGE WITH ENDPOINT INCLUDED arangeWithEndpoint(0,10,2.5) = [0, 2.5, 5, 7.5, 10]
#     returns coordinates in order [theta, r, z]
#
def arangeWithEndpoint(start, stop, step=1, endpoint=True):
	if start == stop:
		return [start]

	arr = np.arange(start, stop, step)
	if endpoint and arr[-1] + step == stop:
		arr = np.concatenate([arr, [stop]])
	return arr



### Setup the simulation
post_proc_only = False

### Switch to configure coarser mesh for faster simulation
coarse_mesh = True

### General parameter setup

# create output directory name based on name of this python script
currDir = os.path.dirname(os.path.abspath(__file__))
modelName = os.path.basename(__file__).replace('.py','')
dataName = modelName + '_output'
Sim_Path = os.path.join(currDir, dataName)
print("Simulation data path: ", Sim_Path)

if not post_proc_only:
    try:
        shutil.rmtree(Sim_Path)           #cleanup
    except:
        pass
    try:
        os.mkdir(Sim_Path)
    except:
        pass


unit = 1e-3 # all length in mm

# 2m band
f0 = 141e6 # simulation center frequency, frequency of interest!
fc = 10e6 # excitation bandwidth for simulation
#fc = 100e6 # excitation bandwidth for simulation
feed_R = 50

plotFrequency = 141e6 # frequency for plotting far field results and return loss
plotBandwidth = 10e6  # bandwidth for return loss plot 

lambda0 = round(C0/f0/unit) # wavelength in mm

### Setup FDTD parameter & excitation function
FDTD = openEMS(EndCriteria=1e-4)
FDTD.SetGaussExcite( f0, fc )
# boundary uses perfectly matched layer in +z direction (beam direction),
# MUR absorbing sheet is good enough for other directions, simulates faster than PML
FDTD.SetBoundaryCond( ['PEC', 'PML_4', 'PML_4', 'PML_4', 'PML_4', 'PML_6'] )  

### Setup Geometry & Mesh
CSX = CSXCAD.ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)

if coarse_mesh:
     cells_per_wavelength = 15
else:
     cells_per_wavelength = 20

max_res = floor(C0 / (f0+fc) / unit / cells_per_wavelength)			# mesh density in air regions etc

# Yagi dimensions
lengths = [ 1038., 995., 959., 949., 935. ]
spaces =  [ 0, 280, 150, 520, 530 ]
space_off = -sum(spaces)/2 # center the antenna in the Z dimension

excitation_length = 5 # mm

# diameter of antenna elements (square)
trace_width = 2

# size of the simulation box
SimBox = np.array([max(lengths)*2, lambda0, sum(spaces)*2])

### Generate properties, primitives and mesh-grid
#initialize the mesh with the "air-box" dimensions
# mesh.AddLine('x', [-SimBox[0]/2, SimBox[0]/2])
mesh.AddLine('x', [ SimBox[0]/2])
mesh.AddLine('y', [-SimBox[1]/2, SimBox[1]/2])
mesh.AddLine('z', [-SimBox[2]/2, SimBox[2]/2])

yagi = CSX.AddMetal('yagi')

for i in range(len(lengths)):
    space_off += spaces[i]  # element offset in z direction

    # port direction is along x-axis
    port_x1 =  0
    port_x2 =  excitation_length

    y1 = - trace_width/2
    y2 =   trace_width/2

    z1 = space_off - trace_width/2
    z2 = space_off + trace_width/2

    if i == 1: # driver
		
        # element_start = [-lengths[i]/2, y1, z1]
        # element_stop  = [port_x1, y2, z2]
        # yagi.AddBox(priority=10,start=element_start,stop=element_stop)

        element_start = [ lengths[i]/2, y1, z1]
        element_stop  = [port_x2, y2, z2]
        yagi.AddBox(priority=10,start=element_start,stop=element_stop)

        port_start = [port_x1,y1,z1]
        port_stop =  [port_x2,y2,z2]
        #port = FDTD.AddLumpedPort(1 ,feed_R, port_start, port_stop, 'x', 1.0, priority=10000)
        port = FDTD.AddLumpedPort(1 ,feed_R/2, port_start, port_stop, 'x', 1.0, priority=10000)

		
    else:
        # element_start = [-lengths[i]/2, y1, z1]
        element_start = [ 0, y1, z1]
        element_stop  = [ lengths[i]/2, y2, z2]
        yagi.AddBox(priority=10,start=element_start,stop=element_stop)
    
    # extra mesh lines near element (but outside metal)
    mesh.AddLine('z', [space_off-2*trace_width, space_off+2*trace_width])

    # extra mesh line near end of elements (but outside metal)
    # mesh.AddLine('x', [-lengths[i]/2-trace_width, -lengths[i]/2+trace_width])
    mesh.AddLine('x', [ lengths[i]/2-trace_width,  lengths[i]/2+trace_width])



# extra mesh lines near element (but outside metal)
mesh.AddLine('y', [-2*trace_width, 2*trace_width])

# extra mesh line in middle of port length
mesh.AddLine('x', [excitation_length/2])


# Add meshing to the "yagi" element
FDTD.AddEdges2Grid(dirs='all', properties=yagi)

# create a smooth mesh between specified fixed mesh lines
if coarse_mesh:
    mesh_stepping_factor = 2
else:
    mesh_stepping_factor = 1.5          

mesh.SmoothMeshLines('all', max_res, ratio=mesh_stepping_factor)

# nf2ff calc
nf2ff_box1 = FDTD.CreateNF2FFBox(opt_resolution=[lambda0/15]*3)

### Run the simulation
if not post_proc_only:
    if 1:  # debugging only
        CSX_file = os.path.join(Sim_Path, modelName + '.xml')
        print("Writing CSXCAD file:" + CSX_file )
        if not os.path.exists(Sim_Path):
            os.mkdir(Sim_Path)
        CSX.Write2XML(CSX_file)
        from CSXCAD import AppCSXCAD_BIN
        os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

    FDTD.Run(Sim_Path, verbose=3, cleanup=False, setup_only=False, debug_pec=False)

### Postprocessing & plotting
numfreq = 501
freq = linspace( plotFrequency-plotBandwidth/2, plotFrequency+plotBandwidth/2, numfreq )
port.CalcPort(Sim_Path, freq,ref_impedance = feed_R)

Zin = port.uf_tot / port.if_tot
s11 = port.uf_ref / port.uf_inc
VSWR = (1+abs(s11))/(1-abs(s11))

# create Touchstone S1P output file using frequency range specified for plot
s1p_name = CSX_file = os.path.join(Sim_Path, modelName + ".s1p")
s1p_file = open(s1p_name, "w")
s1p_file.write("#   Hz   S  RI   R   " + str(feed_R) + "\n")
s1p_file.write("!\n")
for index in range(0, numfreq):
    f = freq[index]
    re = np.real(s11[index])
    im = np.imag(s11[index])
    s1p_file.write(str(f) + " " + str(re) + " " + str(im) + "\n")
s1p_file.close()
print("Wrote Touchstone S1P output file to {}".format(s1p_name))


## Plot the feed point impedance
figure()
plot( freq/1e6, real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' )
grid()
plot( freq/1e6, imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' )
title( 'feed point impedance' )
xlabel( 'frequency (MHz)' )
ylabel( 'impedance ($\Omega$)' )
legend( )


## Plot reflection coefficient S11
figure()
plot( freq/1e6, 20*log10(abs(s11)), 'k-', linewidth=2 )
grid()
title( 'reflection coefficient $S_{11}$' )
xlabel( 'frequency (MHz)' )
ylabel( 'reflection coefficient $|S_{11}|$' )

## Plot VSWR
figure()
plot( freq/1e6, VSWR, 'k-', linewidth=2 )
ylim(1, 10)
grid()
title( 'VSWR' )
xlabel( 'frequency (MHz)' )
ylabel( 'VSWR' )

## Port signals over time

#portVariable.u_data.ui_val[0] = U over time
#portVariable.i_data.ui_val[0] = I over time
#portVariable.u_data.ui_time[0] = timesteps
#portVariable.i_data.ui_time[0] = timesteps (should be same as for voltage)

## Plot port voltage over time
u1 = port.u_data.ui_val[0]
i1 = port.i_data.ui_val[0]
t  = port.u_data.ui_time[0]
figure()
plot(t*1e9,u1, 'k-', label='u1')
xlim(0, 1000)
grid()
title( 'Port voltage vs. time' )
ylabel('Port voltage (V)')
xlabel('Time (ns)')



#######################################################################################################################################
# 2D Farfield plot 
#######################################################################################################################################

# Calculate the far field at phi=0 degrees and at phi=90 degrees
#   Using angles in degrees.
#

print( 'calculating the 2D far field cuts at ' + str(plotFrequency/1E6) + 'MHz')

theta = np.arange(-180.0, 180.0, 2.0)
phi   = [0., 90.]
nf2ff_res = nf2ff_box1.CalcNF2FF(Sim_Path, plotFrequency, theta, phi, center=[0,0,1e-3])


E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])

# rectangular antenna pattern plot 
figure()
plot(theta, np.squeeze(E_norm[:,0]), 'k-', linewidth=2, label='xz-plane')
plot(theta, np.squeeze(E_norm[:,1]), 'r--', linewidth=2, label='yz-plane')
grid()
ylabel('Directivity (dBi)')
xlabel('Theta (deg)')
title('Frequency: {} MHz'.format(plotFrequency/1e6))
legend()

plot_name = os.path.join(Sim_Path, "directivity.png")
savefig(plot_name)

show()

[thliebig]

I currently have no machine to run anything, but looking at your script you should not use PML_4 oder 6. Use Mur if you want quicker results or PML_8 only. Any other than PML_8 is not really tested and I'm not sure how good or if at all they work. If something goes unstable the pml is always the first suspect...

[VolkerMuehlhaus]

That was a great hint, Thorsten, no more convergence issues when using PML_8!

Results look reasonable at first, but I need to check patterns more carefully. There was a message "nf2ff_calc::AddPlane: Error can't determine normal direction..."

Have a great New Years Eve!

[thliebig]

I think you need to disable the xmin direction for the nf2ff recording and calculations.
If everything works as it should, the FF results should be identical...

[VolkerMuehlhaus]

Could you give me an idea how to do that?
I tried the code below, but this is obviously wrong,
TypeError: function() got multiple values for keyword argument 'directions'

nf2ff_box1 = FDTD.CreateNF2FFBox(opt_resolution=[lambda0/15]*3, directions=[False,True,True,True,True,True])

[thliebig]

Oh well, now I confused myself. I had checked 2 days ago that "CreateNF2FFBox" already does it automatically...
But you could check e.g. in AppCSXCAD that no nf2ff planes are created for xmin...
But the nf2ff warning is a bit confusing then.