search

thliebig / openEMS-Project

openEMS is a free and open electromagnetic field solver using the FDTD method.
356 stars 65 forks source link

ISSUE: python AddMSLPort different from octave, not adding Resistor lumped part #124

Closed LubomirJagos closed 8 months ago

LubomirJagos commented 8 months ago

Hi, I'm writing GUI addon for FreeCAD so I have GUI where you define parameters for openEMS simulation and you can generate matlab/octave or python script (now writing python generator) and I was implementing microstrip port and hit this problem:

Question: Octave generated model IS RIGHT AS I EXPECTED AND SHOULD BE, why python IS NOT GENERATING MSL port RESISTANCE?

here is octave code

% OpenEMS FDTD Analysis Automation Script
%
% To be run with GNU Octave or MATLAB.
% FreeCAD to OpenEMS plugin by Lubomir Jagos, 
% see https://github.com/LubomirJagos/FreeCAD-OpenEMS-Export
%
% This file has been automatically generated. Manual changes may be overwritten.
%

close all
clear
clc

%% Change the current folder to the folder of this m-file.
if(~isdeployed)
  mfile_name          = mfilename('fullpath');
  [pathstr,name,ext]  = fileparts(mfile_name);
  cd(pathstr);
end

%% constants
physical_constants;
unit    = 0.001; % Model coordinates and lengths will be specified in mm.
fc_unit = 0.001; % STL files are exported in FreeCAD standard units (mm).

%% switches & options
postprocessing_only = 1;
draw_3d_pattern = 0; % this may take a while...
use_pml = 0;         % use pml boundaries instead of mur

currDir = strrep(pwd(), '\', '\\');
display(currDir);

% --no-simulation : dry run to view geometry, validate settings, no FDTD computations
% --debug-PEC     : generated PEC skeleton (use ParaView to inspect)
openEMS_opts = '--no-simulation';

%% prepare simulation folder
Sim_Path = 'simulation_output';
Sim_CSX = 'Mircostrip_Port_1.xml';
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder

%% setup FDTD parameter & excitation function
max_timesteps = 1000000;
min_decrement = 1e-05; % 10*log10(min_decrement) dB  (i.e. 1E-5 means -50 dB)
FDTD = InitFDTD( 'NrTS', max_timesteps, 'EndCriteria', min_decrement );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BOUNDARY CONDITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BC = {"PEC","PEC","PEC","PEC","PEC","PEC"};
FDTD = SetBoundaryCond( FDTD, BC );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COORDINATE SYSTEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = InitCSX('CoordSystem', 0); % Cartesian coordinate system.
mesh.x = []; % mesh variable initialization (Note: x y z implies type Cartesian).
mesh.y = [];
mesh.z = [];
CSX = DefineRectGrid(CSX, unit, mesh); % First call with empty mesh to set deltaUnit attribute.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXCITATION experimental
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0 = 2.4*1000000000.0;
fc = 1.2*1000000000.0;
FDTD = SetGaussExcite( FDTD, f0, fc );
max_res = c0 / (f0 + fc) / 20;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATERIALS AND GEOMETRY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddMetal( CSX, 'PEC' );

%% MATERIAL - PEC
CSX = AddMetal(CSX, 'PEC');

%% MATERIAL - air
CSX = AddMaterial(CSX, 'air');
CSX = SetMaterialProperty(CSX, 'air', 'Epsilon', 1, 'Mue', 1);

%% MATERIAL - copper sheet
CSX = AddConductingSheet(CSX, 'copper sheet', 5e+07, 35*1e-06);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GRID LINES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% GRID - AIR - air (Fixed Distance)
mesh.x(mesh.x >= -51.0358 & mesh.x <= 56.9278) = [];
mesh.x = [ mesh.x (-51.0358:2:56.9278) ];
mesh.y(mesh.y >= -54.5894 & mesh.y <= 56.6877) = [];
mesh.y = [ mesh.y (-54.5894:2:56.6877) ];
mesh.z(mesh.z >= -40 & mesh.z <= 40) = [];
mesh.z = [ mesh.z (-40:2:40) ];
CSX = DefineRectGrid(CSX, unit, mesh);

%% GRID - XY 20 lines - substrate (Fixed Count)
mesh.x(mesh.x >= -29.5166 & mesh.x <= 30.4928) = [];
mesh.x = [ mesh.x linspace(-29.5166,30.4928,200) ];
mesh.y(mesh.y >= -28.4 & mesh.y <= 30.3971) = [];
mesh.y = [ mesh.y linspace(-28.4,30.3971,200) ];
CSX = DefineRectGrid(CSX, unit, mesh);

%% GRID - Z 5 line - substrate (Fixed Count)
mesh.z(mesh.z >= 0 & mesh.z <= 1.54) = [];
mesh.z = [ mesh.z linspace(0,1.54,5) ];
CSX = DefineRectGrid(CSX, unit, mesh);

%% GRID - Z 5 line - MLS port IN (Fixed Count)
mesh.z(mesh.z >= 0.04 & mesh.z <= 1.54) = [];
mesh.z = [ mesh.z linspace(0.04,1.54,5) ];
CSX = DefineRectGrid(CSX, unit, mesh);

%% GRID - X MSL PORT lines - MLS port IN (Fixed Count)
mesh.x(mesh.x >= 10.6736 & mesh.x <= 14.4867) = [];
mesh.x = [ mesh.x linspace(10.6736,14.4867,200) ];
CSX = DefineRectGrid(CSX, unit, mesh);

%% GRID - Y MSL PORT linex - MLS port IN (Fixed Count)
mesh.y(mesh.y >= -26.8079 & mesh.y <= -3.45277) = [];
mesh.y = [ mesh.y linspace(-26.8079,-3.45277,200) ];
CSX = DefineRectGrid(CSX, unit, mesh);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PORTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
portNamesAndNumbersList = containers.Map();

%% PORT - microstrip sink - MSL port OUT
portStart  = [ 10.9443, 24.184, 1.54 ];
portStop = [ 13.672, -2.33531, 0.04 ];
portR = 48;
portUnits = 1;
portExcitationAmplitude = 3.7;
mslDir = 1;
mslEVec = [0 0 -1]*portExcitationAmplitude;
[CSX port{1}] = AddMSLPort(CSX,9900,1,'copper sheet',portStart, portStop, mslDir, mslEVec, 'FeedShift', 3.0, 'MeasPlaneShift', 4.0, 'Feed_R', portR*portUnits);
portNamesAndNumbersList("MSL port OUT") = 1;

%% PORT - microstrip source - MLS port IN
portStart  = [ 10.6736, -26.8079, 1.54 ];
portStop = [ 14.4867, -3.45277, 0.04 ];
portR = 11;
portUnits = 1;
portExcitationAmplitude = 3.7;
mslDir = 1;
mslEVec = [0 0 -1]*portExcitationAmplitude;
[CSX port{2}] = AddMSLPort(CSX,9800,2,'copper sheet',portStart, portStop, mslDir, mslEVec, 'ExcitePort', true, 'FeedShift', 1.0, 'MeasPlaneShift', 2.0, 'Feed_R', portR*portUnits);
portNamesAndNumbersList("MLS port IN") = 2;

%% PORT - lumped port - lumpedPort_1
portStart = [ -12.8372, 9.2439, -1.7 ];
portStop  = [ 3.54001, 10.6087, 0 ];
portR = 89.0;
portUnits = 1;
portExcitationAmplitude = 2.3;
portDirection = [0 0 1]*portExcitationAmplitude;
[CSX port{3}] = AddLumpedPort(CSX, 10000, 3, portR*portUnits, portStart, portStop, portDirection);
portNamesAndNumbersList("lumpedPort_1") = 3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
CSXGeomPlot( [Sim_Path '/' Sim_CSX] );

if (postprocessing_only==0)
    %% run openEMS
    RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts );
end

here is python code

# OpenEMS FDTD Analysis Automation Script
#
# To be run with python.
# FreeCAD to OpenEMS plugin by Lubomir Jagos, 
# see https://github.com/LubomirJagos/FreeCAD-OpenEMS-Export
#
# This file has been automatically generated. Manual changes may be overwritten.
#
### Import Libraries
import os, tempfile, shutil
from pylab import *
import CSXCAD
from openEMS import openEMS
from openEMS.physical_constants import *
#
#
# Change current path to script file folder
#
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
## constants
unit    = 0.001 # Model coordinates and lengths will be specified in mm.
fc_unit = 0.001 # STL files are exported in FreeCAD standard units (mm).

## switches & options
postprocessing_only = True;
draw_3d_pattern = 0  # this may take a while...
use_pml = 0          # use pml boundaries instead of mur

currDir = os.getcwd()
print(currDir)

# --no-simulation : dry run to view geometry, validate settings, no FDTD computations
# --debug-PEC     : generated PEC skeleton (use ParaView to inspect)
openEMS_opts = '--no-simulation'

## prepare simulation folder
Sim_Path = os.path.join(currDir, 'simulation_output')
Sim_CSX = 'Mircostrip Port 1.xml'
if os.path.exists(Sim_Path):
    shutil.rmtree(Sim_Path)   # clear previous directory
    os.mkdir(Sim_Path)    # create empty simulation folder

## setup FDTD parameter & excitation function
max_timesteps = 1000000
min_decrement = 1e-05 # 10*log10(min_decrement) dB  (i.e. 1E-5 means -50 dB)
CSX = CSXCAD.ContinuousStructure()
FDTD = openEMS(NrTS=max_timesteps, EndCriteria=min_decrement)
FDTD.SetCSX(CSX)

#######################################################################################################################################
# BOUNDARY CONDITIONS
#######################################################################################################################################
BC = ["PEC","PEC","PEC","PEC","PEC","PEC"]
FDTD.SetBoundaryCond(BC)

#######################################################################################################################################
# COORDINATE SYSTEM
#######################################################################################################################################
FDTD.SetCoordSystem(0) # Cartesian coordinate system.
def mesh():
    x,y,z

mesh.x = np.array([]); # mesh variable initialization (Note: x y z implies type Cartesian).
mesh.y = np.array([]);
mesh.z = np.array([]);
openEMS_grid = CSX.GetGrid()
openEMS_grid.SetDeltaUnit(unit) # First call with empty mesh to set deltaUnit attribute.

#######################################################################################################################################
# EXCITATION experimental
#######################################################################################################################################
f0 = 2.4*1000000000.0
fc = 1.2*1000000000.0
FDTD.SetGaussExcite(f0, fc)
max_res = C0 / (f0 + fc) / 20

#######################################################################################################################################
# MATERIALS AND GEOMETRY
#######################################################################################################################################
PEC = CSX.AddMetal('PEC');

## MATERIAL - PEC
materialList = {}
materialList['PEC'] = CSX.AddMetal('PEC')

## MATERIAL - air
materialList = {}
materialList['air'] = CSX.AddMaterial('air')
materialList['air'].SetMaterialProperty(epsilon=1, mue=1)

## MATERIAL - copper sheet
materialList = {}
materialList['copper sheet'] = CSX.AddConductingSheet('copper sheet', conductivity=5e+07, thickness=35*1e-06)

#######################################################################################################################################
# GRID LINES
#######################################################################################################################################

## GRID - AIR - air (Fixed Distance)
mesh.x = np.delete(mesh.x, np.argwhere((mesh.x >= -51.0357) & (mesh.x <= 56.9277)))
mesh.x = np.concatenate((mesh.x, np.arange(-51.0357,56.9277,2)))
mesh.y = np.delete(mesh.y, np.argwhere((mesh.y >= -54.5893) & (mesh.y <= 56.6876)))
mesh.y = np.concatenate((mesh.y, np.arange(-54.5893,56.6876,2)))
mesh.z = np.delete(mesh.z, np.argwhere((mesh.z >= -39.9999) & (mesh.z <= 39.9999)))
mesh.z = np.concatenate((mesh.z, np.arange(-39.9999,39.9999,2)))

## GRID - XY 20 lines - substrate (Fixed Count)
mesh.x = np.delete(mesh.x, np.argwhere((mesh.x >= -29.5165) & (mesh.x <= 30.4927)))
mesh.x = np.concatenate((mesh.x, linspace(-29.5165,30.4927,200)))
mesh.y = np.delete(mesh.y, np.argwhere((mesh.y >= -28.3999) & (mesh.y <= 30.397)))
mesh.y = np.concatenate((mesh.y, linspace(-28.3999,30.397,200)))

## GRID - Z 5 line - substrate (Fixed Count)
mesh.z = np.delete(mesh.z, np.argwhere((mesh.z >= 0.0001) & (mesh.z <= 1.5399)))
mesh.z = np.concatenate((mesh.z, linspace(0.0001,1.5399,5)))

## GRID - Z 5 line - MLS port IN (Fixed Count)
mesh.z = np.delete(mesh.z, np.argwhere((mesh.z >= 0.0401) & (mesh.z <= 1.5399)))
mesh.z = np.concatenate((mesh.z, linspace(0.0401,1.5399,5)))

## GRID - X MSL PORT lines - MLS port IN (Fixed Count)
mesh.x = np.delete(mesh.x, np.argwhere((mesh.x >= 10.6737) & (mesh.x <= 14.4866)))
mesh.x = np.concatenate((mesh.x, linspace(10.6737,14.4866,200)))

## GRID - Y MSL PORT linex - MLS port IN (Fixed Count)
mesh.y = np.delete(mesh.y, np.argwhere((mesh.y >= -26.8078) & (mesh.y <= -3.45288)))
mesh.y = np.concatenate((mesh.y, linspace(-26.8078,-3.45288,200)))

openEMS_grid.AddLine('x', mesh.x)
openEMS_grid.AddLine('y', mesh.y)
openEMS_grid.AddLine('z', mesh.z)

#######################################################################################################################################
# PORTS
#######################################################################################################################################
port = {}
portNamesAndNumbersList = {}
## PORT - microstrip sink - MSL port OUT
portStart  = [ 10.9443, 24.184, 1.54 ]
portStop = [ 13.672, -2.33531, 0.04 ]
portR = 48
portUnits = 1
portExcitationAmplitude = 3.7 * 0
mslDir = 'y'
mslEVec = 'z'
port[1] = FDTD.AddMSLPort(1, materialList['copper sheet'], portStart, portStop, mslDir, mslEVec, excite=portExcitationAmplitude, priority=9900, R=portR*portUnits, FeedShift=3.0, MeasPlaneShift=4.0)
portNamesAndNumbersList["MSL port OUT"] = 1;

## PORT - microstrip source - MLS port IN
portStart  = [ 10.6736, -26.8079, 1.54 ]
portStop = [ 14.4867, -3.45277, 0.04 ]
portR = 11
portUnits = 1
portExcitationAmplitude = 3.7 * 1
mslDir = 'y'
mslEVec = 'z'
port[2] = FDTD.AddMSLPort(2, materialList['copper sheet'], portStart, portStop, mslDir, mslEVec, excite=portExcitationAmplitude, priority=9800, R=portR*portUnits, FeedShift=1.0, MeasPlaneShift=2.0)
portNamesAndNumbersList["MLS port IN"] = 2;

## PORT - lumped port - lumpedPort_1
portStart = [ -12.8372, 9.2439, -1.7 ]
portStop  = [ 3.54001, 10.6087, 0 ]
portR = 89.0
portUnits = 1
portDirection = 'z'
port[3] = FDTD.AddLumpedPort(port_nr=3, R=portR*portUnits, start=portStart, stop=portStop, p_dir=portDirection, priority=10000, excite=0)

#######################################################################################################################################
# RUN
#######################################################################################################################################
### Run the simulation
CSX_file = os.path.join(Sim_Path, Sim_CSX)
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))

if not postprocessing_only:
    FDTD.Run(Sim_Path, verbose=3, cleanup=True)

this is octave generated model which is RIGHT image

here is python generated model which is WRONG image

thliebig commented 8 months ago

Looking at the port definition I think the problem is that for the python interface the keyword is "Feed_R" not just "R".... See here

LubomirJagos commented 8 months ago

yes, that's it, this is missing in python interface documentation for ports here https://docs.openems.de/python/openEMS/ports.html