Merge branch 'master' into multithreading
Sebastian Held
commit
55c0719b4e
|
|
@ -92,7 +92,8 @@ unsigned int Operator::GetNyquistNum(double fmax)
|
|||
{
|
||||
if (dT==0) return 1;
|
||||
double T0 = 1/fmax;
|
||||
return floor(T0/2/dT);
|
||||
m_nyquistTS = floor(T0/2/dT);
|
||||
return m_nyquistTS;
|
||||
}
|
||||
|
||||
bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower)
|
||||
|
|
@ -207,29 +208,35 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
|
|||
double Operator::GetNumberCells()
|
||||
{
|
||||
if (numLines)
|
||||
return (numLines[0]-1)*(numLines[1]-1)*(numLines[2]-1);
|
||||
return (numLines[0])*(numLines[1])*(numLines[2]); //it's more like number of nodes???
|
||||
return 0;
|
||||
}
|
||||
|
||||
void Operator::ShowSize()
|
||||
void Operator::ShowStat()
|
||||
{
|
||||
unsigned int OpSize = 12*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
|
||||
unsigned int FieldSize = 6*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
|
||||
double MBdiff = 1024*1024;
|
||||
|
||||
cout << "---- Stat: FDTD Operator ----" << endl;
|
||||
cout << "Dimensions : " << numLines[0] << "x" << numLines[1] << "x" << numLines[2] << " = " << numLines[0]*numLines[1]*numLines[2] << " Cells (" << numLines[0]*numLines[1]*numLines[2]/1e6 << " MCells)" << endl;
|
||||
cout << "Size of Operator : " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MB) " << endl;
|
||||
cout << "Size of Field-Data: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MB) " << endl;
|
||||
cout << "-----------------------------" << endl;
|
||||
cout << "------- Stat: FDTD Operator -------" << endl;
|
||||
cout << "Dimensions : " << numLines[0] << "x" << numLines[1] << "x" << numLines[2] << " = " << numLines[0]*numLines[1]*numLines[2] << " Cells (" << numLines[0]*numLines[1]*numLines[2]/1e6 << " MCells)" << endl;
|
||||
cout << "Size of Operator : " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MB) " << endl;
|
||||
cout << "Size of Field-Data : " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MB) " << endl;
|
||||
cout << "-----------------------------------" << endl;
|
||||
cout << "Timestep (s) : " << dT << endl;
|
||||
cout << "Nyquist criteria (TS) : " << m_nyquistTS << endl;
|
||||
cout << "Nyquist criteria (s) : " << m_nyquistTS*dT << endl;
|
||||
cout << "Excitation Length (TS) : " << ExciteLength << endl;
|
||||
cout << "Excitation Length (s) : " << ExciteLength*dT << endl;
|
||||
cout << "-----------------------------------" << endl;
|
||||
}
|
||||
|
||||
bool Operator::CalcGaussianPulsExcitation(double f0, double fc)
|
||||
unsigned int Operator::CalcGaussianPulsExcitation(double f0, double fc)
|
||||
{
|
||||
if (dT==0) return false;
|
||||
if (dT==0) return 0;
|
||||
|
||||
ExciteLength = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
|
||||
cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
|
||||
// cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
|
||||
delete[] ExciteSignal;
|
||||
ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
|
||||
ExciteSignal[0]=0.0;
|
||||
|
|
@ -238,16 +245,16 @@ bool Operator::CalcGaussianPulsExcitation(double f0, double fc)
|
|||
ExciteSignal[n] = cos(2.0*PI*f0*(n*dT-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*n*dT/3.0-3,2));
|
||||
// cerr << ExciteSignal[n] << endl;
|
||||
}
|
||||
return true;
|
||||
return GetNyquistNum(f0+fc);
|
||||
}
|
||||
|
||||
bool Operator::CalcSinusExcitation(double f0, int nTS)
|
||||
unsigned int Operator::CalcSinusExcitation(double f0, int nTS)
|
||||
{
|
||||
if (dT==0) return false;
|
||||
if (nTS<=0) return false;
|
||||
if (dT==0) return 0;
|
||||
if (nTS<=0) return 0;
|
||||
|
||||
ExciteLength = (unsigned int)(nTS);
|
||||
cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
|
||||
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
|
||||
delete[] ExciteSignal;
|
||||
ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
|
||||
ExciteSignal[0]=0.0;
|
||||
|
|
@ -256,7 +263,7 @@ bool Operator::CalcSinusExcitation(double f0, int nTS)
|
|||
ExciteSignal[n] = sin(2.0*PI*f0*n*dT);
|
||||
// cerr << ExciteSignal[n] << endl;
|
||||
}
|
||||
return true;
|
||||
return GetNyquistNum(f0);
|
||||
}
|
||||
|
||||
void Operator::DumpOperator2File(string filename)
|
||||
|
|
@ -719,7 +726,7 @@ double Operator::CalcTimestep()
|
|||
}
|
||||
}
|
||||
}
|
||||
cerr << "Operator Timestep: " << dT << endl;
|
||||
// cerr << "Operator Timestep: " << dT << endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -37,10 +37,10 @@ public:
|
|||
|
||||
virtual int CalcECOperator();
|
||||
|
||||
//! Calculate an excitation with center of f0 and the half bandwidth fc
|
||||
virtual bool CalcGaussianPulsExcitation(double f0, double fc);
|
||||
//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps
|
||||
virtual bool CalcSinusExcitation(double f0, int nTS);
|
||||
//! Calculate an excitation with center of f0 and the half bandwidth fc \return number of Nyquist timesteps
|
||||
virtual unsigned int CalcGaussianPulsExcitation(double f0, double fc);
|
||||
//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps \return number of Nyquist timesteps
|
||||
virtual unsigned int CalcSinusExcitation(double f0, int nTS);
|
||||
|
||||
virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
|
||||
virtual void ApplyMagneticBC(bool* dirs);
|
||||
|
|
@ -49,7 +49,7 @@ public:
|
|||
unsigned int GetNyquistNum(double fmax);
|
||||
double GetNumberCells();
|
||||
|
||||
void ShowSize();
|
||||
void ShowStat();
|
||||
|
||||
void DumpOperator2File(string filename);
|
||||
void DumpMaterial2File(string filename);
|
||||
|
|
@ -86,6 +86,7 @@ protected:
|
|||
//Calc timestep only internal use
|
||||
virtual double CalcTimestep();
|
||||
double dT; //FDTD timestep!
|
||||
unsigned int m_nyquistTS;
|
||||
|
||||
//EC elements, internal only!
|
||||
bool Calc_EC();
|
||||
|
|
|
|||
|
|
@ -0,0 +1,28 @@
|
|||
function UI = ReadUI(files, path)
|
||||
|
||||
if (nargin<2)
|
||||
path ='';
|
||||
end
|
||||
|
||||
if (ischar(files))
|
||||
filenames{1}=files;
|
||||
else
|
||||
filenames=files;
|
||||
end
|
||||
|
||||
for n=1:numel(filenames)
|
||||
tmp = load([path filenames{n}]);
|
||||
t = tmp(:,1)';
|
||||
val = tmp(:,2)';
|
||||
|
||||
UI.TD{n}.t = t;
|
||||
UI.TD{n}.val = val;
|
||||
|
||||
dt=t(2)-t(1);
|
||||
val = [val zeros(1,5000)];
|
||||
L=numel(val);
|
||||
UI.FD{n}.f = (0:L-1)/L/dt;
|
||||
UI.FD{n}.f = UI.FD{n}.f(1:floor(L/2));
|
||||
UI.FD{n}.val = fft(val)/L;
|
||||
UI.FD{n}.val = UI.FD{n}.val(1:floor(L/2));
|
||||
end
|
||||
|
|
@ -8,14 +8,18 @@ coil_rad = 10;
|
|||
coil_length = 50;
|
||||
coil_turns = 8;
|
||||
coil_res = 10;
|
||||
port_length = 10;
|
||||
port_length = coil_length/2;
|
||||
port_resist = 50;
|
||||
|
||||
f_max = 50e6;
|
||||
f_excite = 1e9;
|
||||
mesh_size = wire_rad;
|
||||
|
||||
openEMS_Path = [pwd() '/../../']
|
||||
openEMS_opts = '';
|
||||
% openEMS_opts = [openEMS_opts ' --disable-dumps'];
|
||||
% openEMS_opts = [openEMS_opts ' --debug-material'];
|
||||
% openEMS_opts = [openEMS_opts ' --debug-operator'];
|
||||
|
||||
Sim_Path = 'tmp';
|
||||
Sim_CSX = 'helix.xml';
|
||||
|
|
@ -24,15 +28,21 @@ mkdir(Sim_Path);
|
|||
|
||||
%setup FDTD parameter
|
||||
FDTD = InitFDTD(5e5,1e-6);
|
||||
FDTD = SetGaussExcite(FDTD,0.5e9,0.5e9);
|
||||
FDTD = SetGaussExcite(FDTD,f_excite/2,f_excite/2);
|
||||
BC = [1 1 1 1 1 1];
|
||||
FDTD = SetBoundaryCond(FDTD,BC);
|
||||
|
||||
add_Lines = mesh_size * 1.5.^(1:10);
|
||||
add_Lines = add_Lines(find(add_Lines<(3e8/f_excite)/10*1e3));
|
||||
|
||||
%setup CSXCAD geometry
|
||||
CSX = InitCSX();
|
||||
mesh.x = [-35,-25,-20,-15,-12,-11,-10.5,-10,-9.5,-9,-8.5,-8,-7.5,-7,-6.5,-6,-5.5,-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,12,13.5,15,17,18,19,19.5,20,20.5,21,22,23,25,27.5,30,35,45];
|
||||
mesh.y = [-35,-25,-20,-15,-12,-11,-10.5,-10,-9.5,-9,-8.5,-8,-7.5,-7,-6.5,-6,-5.5,-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,12,13,15,17.5,20,25,35];
|
||||
mesh.z = [-25,-15,-10,-5,-2,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5,13,13.5,14,14.5,15,15.5,16,16.5,17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,22.5,23,23.5,24,24.5,25,25.5,26,26.5,27,27.5,28,28.5,29,29.5,30,30.5,31,31.5,32,32.5,33,33.5,34,34.5,35,35.5,36,36.5,37,37.5,38,38.5,39,39.5,40,40.5,41,41.5,42,42.5,43,43.5,44,44.5,45,45.5,46,46.5,47,47.5,48,48.5,49,49.5,50,50.5,51,52,53,55,57.5,60,65,75];
|
||||
mesh.x = -coil_rad-mesh_size : mesh_size : coil_rad+mesh_size+feed_length;
|
||||
mesh.x = [mesh.x(1)-add_Lines mesh.x mesh.x(end)+add_Lines ];
|
||||
mesh.y = -coil_rad-mesh_size : mesh_size : coil_rad+mesh_size;
|
||||
mesh.y = [mesh.y(1)-add_Lines mesh.y mesh.y(end)+add_Lines ];
|
||||
mesh.z = -mesh_size : mesh_size : coil_length+mesh_size;
|
||||
mesh.z = [mesh.z(1)-add_Lines mesh.z mesh.z(end)+add_Lines ];
|
||||
CSX = DefineRectGrid(CSX, 1e-3,mesh);
|
||||
|
||||
%create copper helix and feed lines...
|
||||
|
|
@ -43,7 +53,13 @@ CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6);
|
|||
dt = 1.0/coil_res;
|
||||
height=0;
|
||||
wire.Vertex = {};
|
||||
count=0;
|
||||
p(1,1) = coil_rad + feed_length;
|
||||
p(2,1) = 0;
|
||||
p(3,1) = 0.5*(coil_length-port_length);
|
||||
p(1,2) = coil_rad + feed_length;
|
||||
p(2,2) = 0;
|
||||
p(3,2) = 0;
|
||||
count=2;
|
||||
for n=0:coil_turns-1
|
||||
for m=0:coil_res
|
||||
count = count + 1;
|
||||
|
|
@ -53,22 +69,21 @@ for n=0:coil_turns-1
|
|||
end
|
||||
height = height + coil_length/coil_turns;
|
||||
end
|
||||
p(1,count+1) = coil_rad + feed_length;
|
||||
p(2,count+1) = 0;
|
||||
p(3,count+1) = coil_length;
|
||||
p(1,count+2) = coil_rad + feed_length;
|
||||
p(2,count+2) = 0;
|
||||
p(3,count+2) = 0.5*(coil_length+port_length);
|
||||
CSX = AddWire(CSX, 'copper', 0, p, wire_rad);
|
||||
|
||||
start = [coil_rad, 0 , 0];stop = [coil_rad+feed_length, 0 , 0];
|
||||
CSX = AddCylinder(CSX,'copper',0 ,start,stop,wire_rad);
|
||||
|
||||
start(3)=coil_length;stop(3)=coil_length;
|
||||
CSX = AddCylinder(CSX,'copper',0 ,start,stop,wire_rad);
|
||||
|
||||
start(3)=0;start(1)=coil_rad+feed_length;
|
||||
CSX = AddCylinder(CSX,'copper',0 ,start,stop,wire_rad);
|
||||
|
||||
CSX = AddMaterial(CSX,'resist');
|
||||
kappa = port_length/port_resist/wire_rad^2/pi/1e-3
|
||||
kappa = port_length/port_resist/wire_rad^2/pi/1e-3;
|
||||
CSX = SetMaterialProperty(CSX,'resist','Kappa',kappa);
|
||||
|
||||
start(3)=(coil_length-port_length)/2;stop(3)=(coil_length+port_length)/2;
|
||||
start=[coil_rad+feed_length 0 (coil_length-port_length)/2];
|
||||
stop=[coil_rad+feed_length 0 (coil_length+port_length)/2];
|
||||
%start(3)=(coil_length-port_length)/2;stop(3)=(coil_length+port_length)/2;
|
||||
CSX = AddCylinder(CSX,'resist',5 ,start,stop,wire_rad);
|
||||
|
||||
%excitation
|
||||
|
|
@ -77,7 +92,7 @@ CSX = AddCylinder(CSX,'excite', 0 ,start,stop,wire_rad);
|
|||
|
||||
%voltage calc
|
||||
CSX = AddProbe(CSX,'ut1',0);
|
||||
CSX = AddBox(CSX,'ut1', 0 ,start,stop);
|
||||
CSX = AddBox(CSX,'ut1', 0 ,stop,start);
|
||||
|
||||
%current calc
|
||||
CSX = AddProbe(CSX,'it1',1);
|
||||
|
|
@ -86,6 +101,17 @@ start(1) = start(1)-2;start(2) = start(2)-2;
|
|||
stop(1) = stop(1)+2;stop(2) = stop(2)+2;
|
||||
CSX = AddBox(CSX,'it1', 0 ,start,stop);
|
||||
|
||||
%dump
|
||||
CSX = AddDump(CSX,'Et_',0,0);
|
||||
start = [mesh.x(1) , 0 , mesh.z(1)];
|
||||
stop = [mesh.x(end) , 0 , mesh.z(end)];
|
||||
CSX = AddBox(CSX,'Et_',0 , start,stop);
|
||||
|
||||
CSX = AddDump(CSX,'Ht_',1,0);
|
||||
start = [mesh.x(1) , 0 , mesh.z(1)];
|
||||
stop = [mesh.x(end) , 0 , mesh.z(end)];
|
||||
CSX = AddBox(CSX,'Ht_',0 , start,stop);
|
||||
|
||||
%Write openEMS compatoble xml-file
|
||||
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
|
||||
|
||||
|
|
@ -97,3 +123,26 @@ disp(command);
|
|||
system(command)
|
||||
cd(savePath);
|
||||
|
||||
%%%post-proc
|
||||
U = ReadUI('ut1','tmp/');
|
||||
I = ReadUI('it1','tmp/');
|
||||
|
||||
Z = U.FD{1}.val./I.FD{1}.val;
|
||||
f = U.FD{1}.f;
|
||||
L = imag(Z)./(f*2*pi);
|
||||
R = real(Z);
|
||||
ind = find(f<f_max);
|
||||
|
||||
subplot(2,1,1);
|
||||
plot(f(ind)*1e-6,L(ind)*1e9,'Linewidth',2);
|
||||
xlabel('frequency (MHz)');
|
||||
ylabel('coil inductance (nH)');
|
||||
grid on;
|
||||
subplot(2,1,2);
|
||||
plot(f(ind)*1e-6,R(ind),'Linewidth',2);
|
||||
xlabel('frequency (MHz)');
|
||||
ylabel('resistance (\Omega)');
|
||||
grid on;
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,11 @@
|
|||
#!/bin/bash
|
||||
|
||||
#clear LD_LIBRARY_PATH
|
||||
export LD_LIBRARY_PATH=
|
||||
|
||||
exec /home/thorsten/devel/openEMS/openEMS $@
|
||||
#get path to openEMS
|
||||
openEMS_PATH=`dirname $0`
|
||||
|
||||
#execute openEMS
|
||||
exec $openEMS_PATH/openEMS $@
|
||||
|
||||
|
|
|
|||
15
openems.cpp
15
openems.cpp
|
|
@ -198,10 +198,11 @@ int openEMS::SetupFDTD(const char* file)
|
|||
if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1);
|
||||
|
||||
FDTD_Op->CalcECOperator();
|
||||
|
||||
unsigned int Nyquist = 0;
|
||||
if (Excit_Type==0)
|
||||
{
|
||||
if (!FDTD_Op->CalcGaussianPulsExcitation(f0,fc))
|
||||
Nyquist = FDTD_Op->CalcGaussianPulsExcitation(f0,fc);
|
||||
if (!Nyquist)
|
||||
{
|
||||
cerr << "openEMS: excitation setup failed!!" << endl;
|
||||
exit(2);
|
||||
|
|
@ -209,7 +210,8 @@ int openEMS::SetupFDTD(const char* file)
|
|||
}
|
||||
else if (Excit_Type==1)
|
||||
{
|
||||
if (!FDTD_Op->CalcSinusExcitation(f0,NrTS))
|
||||
Nyquist = FDTD_Op->CalcSinusExcitation(f0,NrTS);
|
||||
if (!Nyquist)
|
||||
{
|
||||
cerr << "openEMS: excitation setup failed!!" << endl;
|
||||
exit(2);
|
||||
|
|
@ -232,13 +234,10 @@ int openEMS::SetupFDTD(const char* file)
|
|||
|
||||
time_t OpDoneTime=time(NULL);
|
||||
|
||||
FDTD_Op->ShowSize();
|
||||
FDTD_Op->ShowStat();
|
||||
|
||||
FDTD_Op->ApplyMagneticBC(PMC);
|
||||
|
||||
cout << "Nyquist number of timesteps: " << FDTD_Op->GetNyquistNum(f0+fc) << endl;
|
||||
unsigned int Nyquist = FDTD_Op->GetNyquistNum(f0+fc);
|
||||
|
||||
cout << "Creation time for operator: " << difftime(OpDoneTime,startTime) << " s" << endl;
|
||||
|
||||
//create FDTD engine
|
||||
|
|
@ -331,7 +330,7 @@ int openEMS::SetupFDTD(const char* file)
|
|||
|
||||
void openEMS::RunFDTD()
|
||||
{
|
||||
cout << "Running FDTD engine... this may take a while... grab a coup of coffee?!?" << endl;
|
||||
cout << "Running FDTD engine... this may take a while... grab a cup of coffee?!?" << endl;
|
||||
|
||||
ProcessFields ProcField(FDTD_Op,FDTD_Eng);
|
||||
double maxE=0,currE=0;
|
||||
|
|
|
|||
Loading…
Reference in New Issue