Source driven problem using Method of Moments

Applications
1

Hello,

I am attempting to solve the following problem where I have a microstrip, and I excite the edge m using method of moments (MoM) . This is related to a previous post but it is a different question, which is why I am making a new thread.
Microstrip
I create a set of zero coefficients and place a voltage at the DOF point of my edge. Using these coefficients, I define a Gridfunction and solve the system. Where I am solving for the Electric Field Integral Equation (EFIE), represented by the maxwell single layer operator.

Once solved, I try to find the input impedance of the microstrip. Following the delta-gap excitation definition, we can relate the voltage to the edge’s basis function, the edge length, and current across the edge. Because the normal component of the RWG function is always one, then we obtain Z_{in} = V_{in} / l_m^2 I_m. Giving us the following result for the impedance.

newMicrostrip
newMicrostrip1000×800 55.9 KB

Compared to the literature (from W.C.Gibson MoM in Electromagnetics book):
StripImp
The magnitude is off, and the imaginary component is flipped. Something is wrong, but I do not know what. I need help identifying the problem.

Thank you!

Here is the code for the problem:

import bempp
import numpy as np
from bempp.api.linalg import lu
from scipy.constants.constants import speed_of_light

#Set-up problem
stepSize = 55; #Solution Step
fMax = 1000E6 #1Ghz
fMin = 100E6  #100Mhz 
fStep = (fMax - fMin) / (stepSize-1)
frequency = np.arange(stepSize)*fStep + fMin
vacuum_permittivity = 8.854187817E-12
vacuum_permeability = 4 * np.pi * 1E-7
w = 2 * np.pi * frequency
k = w * np.sqrt(vacuum_permittivity * vacuum_permeability)
wavelength = speed_of_light / fMax
Vin = 1.0 #Excitation

#Rectangle Microstrip
W = 0.003 # 3mm
L = 1 #1m
strip = np.array([[-L/2, W/2, 0], [L/2, W/2, 0], [L/2, -W/2, 0], [-L/2, -W/2, 0]])

#Create empty Zin
Zin = np.zeros(stepSize, dtype = 'complex')

#Mesh
grid = bempp.api.shapes.shapes.screen(strip, h=(1/10)*wavelength)

edgeInd = 70 #Edge of interest          

#Basis Functions to represent problem
div_space = bempp.api.function_space(grid, "RWG", 0)
curl_space = bempp.api.function_space(grid, "SNC", 0)

#Set-Up Coefficients
coefficients=np.zeros(div_space.global_dof_count, dtype = 'complex')
for edges, dofs, mults in zip(grid.element_edges.T, div_space.local2global, div_space.local_multipliers): 
    for i, e in enumerate(edges):
        if e == edgeInd:
            place = dofs[i]
coefficients[place] = Vin #Place excitation

#Excitation GridFunction
trace_fun = bempp.api.GridFunction(div_space, coefficients= coefficients, dual_space=curl_space)
   
for ite in range(stepSize):
    #Set-up Appropriate Operators - EFIE
    efie = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, k[ite])
    
    #Solve System
    solved_system = lu(efie, trace_fun)
    
    #Impedance
    I = solved_system.coefficients[place]*W*W
    Zin[ite] = Vin / I
   
#Plot   
import matplotlib
from matplotlib import pylab as plt
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0) 

plt.plot(frequency*1E-9, np.real(Zin), label = 'Resistance')
plt.plot(frequency*1E-9, np.imag(Zin), '--', label = 'Reactance')
plt.xlim((np.min(frequency)*1E-9,np.max(frequency)*1E-9))
plt.xlabel('Frequency')
plt.ylabel('Impedance')
plt.title('Impedance of Microstrip')
plt.grid(linewidth=1)
plt.legend()
plt.show()
2

Hi. Without further information I won’t be able to help. We are happy to to give advise on the software. But I am not familiar enough with your application to give competent advise on what you are doing there.

3

Okay, could you help me with computing the surface current at a certain point(s)/ or edge of my grid once I have solved the system?

For example, in this case, I am interested in obtaining the current, I_m, across the edge m. Because with this current I can solve for the impedance, using the equation I posted earlier, Z_{in} = V_{in} / l_m^2 I_m.

4

Hello,

What operator does BEMPP have that can evaluate the current of my grid?

For example, in the problem where I have \boldsymbol{F} = \int f \cdot u, where f is the mathematical description of my geometry, u is the complex scalar values of my solution data, and \boldsymbol{F} is the quantity I am solving for, in my case current.

5

Hello,

I tried to reproduce this results without exact success. Using the same code above I get the following results. They look inverted and with some glitches:

image
image868×701 57.3 KB

6

I achieved a more coherent results using the following code:

#Set-up problem
stepSize = 55 #Solution Step
fMax = 1000e6 #1Ghz
fMin = 100e6  #100Mhz 
fStep = (fMax - fMin) / (stepSize-1)
frequency = np.arange(stepSize)*fStep + fMin
vacuum_permittivity = 8.854187817e-12
vacuum_permeability = 4 * np.pi * 1e-7
w = 2 * np.pi * frequency
k = w * np.sqrt(vacuum_permittivity * vacuum_permeability)
wavelength = speed_of_light / fMax
Vin = 1.0 #Excitation

#Rectangle Microstrip
W = 0.003 # 3mm
L = 1.0 #1m
strip = np.array([[-L/2.0, W/2.0, 0.0], [L/2.0, W/2.0, 0.0], [L/2.0, -W/2.0, 0.0], [-L/2.0, -W/2.0, 0.0]])

#Create empty Zin
Zin = np.zeros(stepSize, dtype = 'complex')

#Mesh
GridSizeh=0.1*wavelength
grid = bempp.api.shapes.shapes.screen(strip, h=GridSizeh)

edgeInd = 70 #Edge of interest          

#Basis Functions to represent problem
div_space = bempp.api.function_space(grid, "RWG", 0)
curl_space = bempp.api.function_space(grid, "SNC", 0)

#Set-Up Coefficients
coeff=np.zeros(div_space.global_dof_count, dtype = 'complex')
"""
#for edges, dofs, mults in zip(grid.element_edges.T, div_space.local2global, div_space.local_multipliers):
    for i, e in enumerate(edges):
        if e == edgeInd:
            place = dofs[i]
"""
pointV1 = [0.0, W/2.0, 0.0]
place=BemppIdentifyDOFcellsFromPoint(pointV1,grid,div_space,GridSizeh,False)

coeff[place] = Vin #Place excitation

#Excitation GridFunction
trace_fun = bempp.api.GridFunction(div_space, coefficients= coeff, dual_space=curl_space)
#bempp.api.enable_console_logging()
for ite in range(stepSize):    
    #Set-up Appropriate Operators - EFIE
    efie = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, k[ite])#, assembler="dense",device_interface="opencl")
    
    #Solve System
    #lu_factorOperator=bempp.api.linalg.direct_solvers.compute_lu_factors(efie)
    solved_system = bempp.api.linalg.lu(efie, trace_fun)#, lu_factor=lu_factorOperator)
    
    #Impedance
    I = solved_system.coefficients[place]*W*W
    Zin[ite] = Vin / I
   
#Plot   
import matplotlib
from matplotlib import pylab as plt
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0) 

plt.plot(frequency*1e-9, np.real(Zin), label = 'Resistance')
plt.plot(frequency*1e-9, np.imag(Zin), '--', label = 'Reactance')
plt.xlim((np.min(frequency)*1e-9,np.max(frequency)*1e-9))
plt.xlabel('Frequency')
plt.ylabel('Impedance')
plt.title('Impedance of Microstrip')
plt.grid(linewidth=1)
plt.legend()
7

where:

def BemppIdentifyDOFcellsFromPoint(pointV1aux,gridAux,div_spaceAux,GridSizehAux,CelldofAux=False):
  ## Seems better for coefficient excitation to excite a single edge, and correspondingly its DOF. Variable CelldofAux=False identifies the closests edge DOF. Variable CelldofAux=True identifies the cells DOFs
  #Get edges
  EdgeInd=[] # Number identifying the edge number. A list because maybe more than one point is found
  v1EdgeInd=[] # List of vertice 1
  v2EdgeInd=[] # List of vertice 2
  Diff1EdgeInd=[] # List of distances of vertice 1
  Diff2EdgeInd=[] # List of distances of vertice 2
  for i in np.arange(0,gridAux.number_of_edges,1): # First identify a few candidate edges
    v1=gridAux.edges[0,i] # vertice 1 of the edge
    V1axis=np.abs(gridAux.vertices[:, v1]-pointV1aux)
    DiffPoint1=np.sqrt(V1axis[0]**2+V1axis[1]**2+V1axis[2]**2)
    v2=gridAux.edges[1,i] # vertice 2 of the edge
    V2axis=np.abs(gridAux.vertices[:, v2]-pointV1aux)
    DiffPoint2=np.sqrt(V2axis[0]**2+V2axis[1]**2+V2axis[2]**2)
    #if (DiffPoint1<1.5*GridSizehAux and DiffPoint2<1.5*GridSizehAux):
    EdgeInd.append(i)
    Diff1EdgeInd.append(DiffPoint1)
    Diff2EdgeInd.append(DiffPoint2)
    v1EdgeInd.append(gridAux.vertices[:, v1])
    v2EdgeInd.append(gridAux.vertices[:, v2])
  # Order from minimum distance (the sum of the two distances between the two vertices of the edge) to maximum distance the identified edges
  EdgeInd=np.array(EdgeInd)# Convert to numpy array
  v1EdgeInd=np.array(v1EdgeInd)# Convert to numpy array
  v2EdgeInd=np.array(v2EdgeInd)# Convert to numpy array
  Diff1EdgeInd=np.array(Diff1EdgeInd)# Convert to numpy array
  Diff2EdgeInd=np.array(Diff2EdgeInd)# Convert to numpy array
  TotalDiffEdgeInd=Diff1EdgeInd+Diff2EdgeInd # Total distance of the edge vertices to the point
  OrderIndicesEdgeInd=np.argsort(TotalDiffEdgeInd)
  OrderedMinEdgeIndPot=EdgeInd[OrderIndicesEdgeInd]
  OrderedMinv1EdgeIndPot=v1EdgeInd[OrderIndicesEdgeInd]
  OrderedMinv2EdgeIndPot=v2EdgeInd[OrderIndicesEdgeInd]
  OrderedMinDiff1EdgeIndPot=Diff1EdgeInd[OrderIndicesEdgeInd]
  OrderedMinDiff2EdgeIndPot=Diff2EdgeInd[OrderIndicesEdgeInd]
  # Consider the first two edges as the rulers. Now, order the other edges that share vertices
  # Find the common vertices of the first two edges, and then find the opposite vertices that the third edge should have
  OrderedMinEdgeInd=np.zeros_like(OrderedMinEdgeIndPot)
  OrderedMinv1EdgeInd=np.zeros_like(OrderedMinv1EdgeIndPot)
  OrderedMinv2EdgeInd=np.zeros_like(OrderedMinv2EdgeIndPot)
  OrderedMinDiff1EdgeInd=np.zeros_like(OrderedMinDiff1EdgeIndPot)
  OrderedMinDiff2EdgeInd=np.zeros_like(OrderedMinDiff2EdgeIndPot)
  if (CelldofAux):
    OrderedMinEdgeInd[0:2]=OrderedMinEdgeIndPot[0:2]
    OrderedMinv1EdgeInd[:,0:2]=OrderedMinv1EdgeIndPot[:,0:2]
    OrderedMinv2EdgeInd[:,0:2]=OrderedMinv2EdgeIndPot[:,0:2]
    OrderedMinDiff1EdgeInd[0:2]=OrderedMinDiff1EdgeIndPot[0:2]
    OrderedMinDiff2EdgeInd[0:2]=OrderedMinDiff2EdgeIndPot[0:2]

    v1a=gridAux.edges[0,OrderedMinEdgeInd[0]] # vertice a of edge 1
    v1b=gridAux.edges[1,OrderedMinEdgeInd[0]] # vertice b of edge 1
    v2a=gridAux.edges[0,OrderedMinEdgeInd[1]] # vertice a of edge 2
    v2b=gridAux.edges[1,OrderedMinEdgeInd[1]] # vertice b of edge 2
    if (v1a==v2a):
      v3a=v1b
      v3b=v2b
    elif(v1a==v2b):
      v3a=v1b
      v3b=v2a
    elif(v1b==v2a):
      v3a=v1a
      v3b=v2b
    else: # v1b==v2b
      v3a=v1a
      v3b=v2a
    # Find the third edge
    Edge3Ind=0
    for iIter in range(2,len(OrderedMinEdgeInd),1):
      if (((v3a==gridAux.edges[0,OrderedMinEdgeIndPot[iIter]]).all() and (v3b==gridAux.edges[1,OrderedMinEdgeIndPot[iIter]]).all() ) or ((v3b==gridAux.edges[0,OrderedMinEdgeIndPot[iIter]]).all()  and (v3a==gridAux.edges[1,OrderedMinEdgeIndPot[iIter]]).all() )):
        Edge3Ind=iIter
    OrderedMinEdgeInd[2]=OrderedMinEdgeIndPot[Edge3Ind]
    OrderedMinv1EdgeInd[2,:]=OrderedMinv1EdgeIndPot[Edge3Ind,:]
    OrderedMinv2EdgeInd[2,:]=OrderedMinv2EdgeIndPot[Edge3Ind,:]
    OrderedMinDiff1EdgeInd[2]=OrderedMinDiff1EdgeIndPot[Edge3Ind]
    OrderedMinDiff2EdgeInd[2]=OrderedMinDiff2EdgeIndPot[Edge3Ind]
    # Find a fourth edge, following the min distance criterion
    if (Edge3Ind==2):
      OrderedMinEdgeInd[3]=OrderedMinEdgeIndPot[3]
      OrderedMinv1EdgeInd[3,:]=OrderedMinv1EdgeIndPot[3,:]
      OrderedMinv2EdgeInd[3,:]=OrderedMinv2EdgeIndPot[3,:]
      OrderedMinDiff1EdgeInd[3]=OrderedMinDiff1EdgeIndPot[3]
      OrderedMinDiff2EdgeInd[3]=OrderedMinDiff2EdgeIndPot[3]
    else:
      OrderedMinEdgeInd[3]=OrderedMinEdgeIndPot[2]
      OrderedMinv1EdgeInd[3,:]=OrderedMinv1EdgeIndPot[2,:]
      OrderedMinv2EdgeInd[3,:]=OrderedMinv2EdgeIndPot[2,:]
      OrderedMinDiff1EdgeInd[3]=OrderedMinDiff1EdgeIndPot[2]
      OrderedMinDiff2EdgeInd[3]=OrderedMinDiff2EdgeIndPot[2]
    #print('OrderedMinEdgeInd: '+str(OrderedMinEdgeInd))
    # Should take the first three edges, that normally define a cell. In can happen that the point of inteest falls in and edge or vertice (then this situation should be handled)
    # But take the first six edges, because usually at the grid border less dofs are identified
    NpotentialEdges=4
    EdgeIndCell=OrderedMinEdgeInd[0:NpotentialEdges]
    v1IndCell=OrderedMinv1EdgeInd[0:NpotentialEdges]
    v2IndCell=OrderedMinv2EdgeInd[0:NpotentialEdges]
    EdgeDiff1IndCell=OrderedMinDiff1EdgeInd[0:NpotentialEdges]
    EdgeDiff2IndCell=OrderedMinDiff2EdgeInd[0:NpotentialEdges]
    # Checks
    if (np.sum(EdgeDiff1IndCell==0.0)>0 or np.sum(EdgeDiff2IndCell==0.0)>0): # Hit vertice
      print("Vertice hit. It should be dealt a side!")
    DotProdArray1=np.zeros((3),dtype=np.float32)
    DotProdArray2=np.zeros((3),dtype=np.float32)
    AbsProdArray1=np.zeros((3),dtype=np.float32)
    AbsProdArray2=np.zeros((3),dtype=np.float32)
    for iIter in range(0,3,1):
      DotProdArray1[iIter]=np.abs(np.dot(pointV1aux-v1IndCell[iIter],v2IndCell[iIter]-v1IndCell[iIter]))
      DotProdArray2[iIter]=np.abs(np.dot(pointV1aux-v2IndCell[iIter],v2IndCell[iIter]-v1IndCell[iIter]))
      AbsProdArray1[iIter]=np.sqrt((pointV1aux[0]-v1IndCell[iIter,0])**2+(pointV1aux[1]-v1IndCell[iIter,1])**2+(pointV1aux[2]-v1IndCell[iIter,2])**2)*np.sqrt((v2IndCell[iIter,0]-v1IndCell[iIter,0])**2+(v2IndCell[iIter,1]-v1IndCell[iIter,1])**2+(v2IndCell[iIter,2]-v1IndCell[iIter,2])**2)
      AbsProdArray2[iIter]=np.sqrt((pointV1aux[0]-v2IndCell[iIter,0])**2+(pointV1aux[1]-v2IndCell[iIter,1])**2+(pointV1aux[2]-v2IndCell[iIter,2])**2)*np.sqrt((v2IndCell[iIter,0]-v1IndCell[iIter,0])**2+(v2IndCell[iIter,1]-v1IndCell[iIter,1])**2+(v2IndCell[iIter,2]-v1IndCell[iIter,2])**2)
    # Probably now hitting an edge is well dealt
    #if (np.sum(DotProdArray1==AbsProdArray1)>0 or np.sum(DotProdArray2==AbsProdArray2)>0): # Hit edge
    #  print("Edge hit. It should be dealt a side!")
    # Get DOF
  else:
    EdgeIndCell=OrderedMinEdgeIndPot[0:1]
    v1IndCell=OrderedMinv1EdgeIndPot[0:1]
    v2IndCell=OrderedMinv2EdgeIndPot[0:1]
    EdgeDiff1IndCell=OrderedMinDiff1EdgeIndPot[0:1]
    EdgeDiff2IndCell=OrderedMinDiff2EdgeIndPot[0:1]
  
  Potentialdof=[]
  """
  dofAux=pointInd
  for iIter in range(0,len(EdgeIndCell),1):
    for edges, dofs in zip(gridAux.edges, div_spaceAux.local2global): 
      for e, d in zip(edges, dofs): 
        if e == EdgeIndCell[iIter]: 
          dof[iIter] = d
  """
  #Set-Up Coefficients
  for iIter in range(0,len(EdgeIndCell),1):
    for edges, dofs in zip(gridAux.element_edges.T, div_spaceAux.local2global):#, div_spaceAux.local_multipliers): 
      for i, e in enumerate(edges):
        if (e == EdgeIndCell[iIter]):
          Potentialdof.append(dofs[i])
  if (CelldofAux): # DOFs forming a cell
    dofAux = ([*set(Potentialdof)])[0:3] # Remove duplicate DOFs from the list, and consider the first three dofs (to form a cell)
  else: # DOF of the closest edge
    dofAux = [([*set(Potentialdof)])[0]] # Remove duplicate DOFs from the list, and consider only the first dofs (of an edge)
  #print('dof: '+str(dofAux))
  return dofAux

image
image868×701 54.3 KB